{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy.optimize import fsolve\n",
    "from scipy.optimize import root\n",
    "from scipy.optimize import least_squares\n",
    "import math\n",
    "from scipy.misc import derivative\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import rc\n",
    "rc('text',usetex=True)\n",
    "from pandas import DataFrame\n",
    "import matplotlib.cm as cm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Online Appendix -- Learning about Usability "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "tolerance = 1e-8 \n",
    "\n",
    "class Usability:\n",
    "    def __init__(self,\n",
    "                 given_beta = 0.996,      #discount factor \n",
    "                 given_alpha_hat = 0.005, \n",
    "                 given_alpha_u = 0.99, #0.77\n",
    "                 given_alpha_l = 0.01,\n",
    "                 given_c = 0.025, \n",
    "                 given_s = 0.58,        #chance of having alpha_cH when state is H\n",
    "                 given_eta = 0.7,        #prob of getting signals \n",
    "                 given_chi = 0.07,      #prob of getting public signals \n",
    "                 given_B = 0.97      #DM utility function parameter \n",
    "                ):\n",
    "        self.beta = given_beta\n",
    "        self.alpha_hat = given_alpha_hat\n",
    "        self.alpha_u = given_alpha_u\n",
    "        self.alpha_l = given_alpha_l\n",
    "        self.c = given_c\n",
    "        self.s = given_s \n",
    "        self.eta = given_eta \n",
    "        self.chi = given_chi \n",
    "        self.B = given_B \n",
    "\n",
    "    def equm_allocation(self):\n",
    "        beta = self.beta\n",
    "        alpha_hat = self.alpha_hat \n",
    "        alpha_u = self.alpha_u\n",
    "        alpha_l = self.alpha_l\n",
    "        c = self.c\n",
    "        s = self.s\n",
    "        eta = self.eta\n",
    "        chi = self.chi\n",
    "        B = self.B\n",
    "        \n",
    "        zeta = eta - chi\n",
    "        ell_lowerbar = 1e-9\n",
    "                \n",
    "        def utility(q): \n",
    "            return B*math.log(max(q,ell_lowerbar)) - B*math.log(ell_lowerbar) + q\n",
    "        \n",
    "        def S(ell): # q = ell \n",
    "            return B*math.log(max(ell,ell_lowerbar)) - B*math.log(ell_lowerbar)\n",
    "        \n",
    "        def S_prime(ell):\n",
    "            return B/ell\n",
    "\n",
    "        q_star = math.inf \n",
    "        w_star = q_star # when buyer has all the bargaining power\n",
    "\n",
    "        # closed form solution for pi_d\n",
    "        alpha_H = (1-chi)*alpha_hat + chi*(s*alpha_u + (1-s)*alpha_l)\n",
    "        alpha_L = (1-chi)*alpha_hat + chi*(s*alpha_l + (1-s)*alpha_u)\n",
    "        \n",
    "        if alpha_L>alpha_H:\n",
    "            print(\"Error: Wrong alpha_H and alpha_L, alpha_L>alpha_H\") \n",
    "            \n",
    "        if B*alpha_H - c < 0:\n",
    "            print(\"Error: Wrong alpha_H\", alpha_H) \n",
    "            \n",
    "        #if B*alpha_L - c > 0:\n",
    "        #    print(\"Error: wrong alpha_L\",alpha_L)\n",
    "        \n",
    "        Gamma = ( (1-beta*(1-eta))/(2*s*beta*eta)\n",
    "                - math.sqrt(((1-beta*(1-eta))/(2*s*beta*eta))**2 - (1-s)/s))\n",
    "\n",
    "        k = (math.log(Gamma) / math.log((1-s)/s))\n",
    "        print('k',k)\n",
    "        \n",
    "        if k<= 0: \n",
    "            print('Error: k<=0')\n",
    "            \n",
    "        pi_d_cf = ((((1-beta)/(beta*eta)+s*(1-Gamma))*c-s*B*alpha_L*(1-Gamma)-(c-alpha_L*B)*(2*s-1))\n",
    "            /(s*B*(alpha_H-alpha_L)*(1-Gamma)-(c-alpha_L*B)*(2*s-1)))\n",
    "        \n",
    "        if pi_d_cf >=1:\n",
    "            print(\"Error: pi_d is larger than 1\")\n",
    "            \n",
    "        if pi_d_cf <0:\n",
    "            print(\"Error: pi_d is negative\")\n",
    "        \n",
    "        def update(y): \n",
    "            \" belief updating after receiving a good news, get pi_prime given pi \"\n",
    "            return y*s/(y*s + (1-y)*(1-s))       \n",
    "        \n",
    "        diff = 10 #to calculate the time horizon\n",
    "        pi = pi_d_cf \n",
    "        count = 2\n",
    "        while diff > tolerance:\n",
    "            pi_prime = update(pi)\n",
    "            count += 1 \n",
    "            diff = 1 - pi_prime\n",
    "            pi = pi_prime\n",
    "        n = count #time horizon grid number\n",
    "        \n",
    "        a = np.zeros(n)\n",
    "        pi_array = np.copy(a) \n",
    "        pi_array[1] = pi_d_cf\n",
    "        sbar_array = np.copy(a) \n",
    "        alpha_array = np.copy(a) \n",
    "            \n",
    "        pi = np.copy(pi_d_cf)  \n",
    "        sbar_array[1] = pi_d_cf*s + (1-pi_d_cf)*(1-s) \n",
    "        sbar_array[0] = 1-s\n",
    "        alpha_array[0] = alpha_L\n",
    "        alpha_array[1] = pi_d_cf*alpha_H + (1-pi_d_cf)*alpha_L\n",
    "        \n",
    "        for j in range(1,n-1):\n",
    "            pi_prime = update(pi) \n",
    "            pi_array[j+1] = pi_prime \n",
    "            sbar_array[j+1] = pi_prime*s + (1-pi_prime)*(1-s)  \n",
    "            alpha_array[j+1] = pi_prime*alpha_H + (1-pi_prime)*alpha_L\n",
    "            pi = pi_prime   \n",
    "        \n",
    "        phi_h = (beta*B*alpha_H - beta*c)/(1-beta)\n",
    "            \n",
    "        phi_l = (beta*B*alpha_L - beta*c)/(1-beta)\n",
    "                \n",
    "        phi_hat_array = beta*(B*alpha_array - c)/(1-beta) # stationary price \n",
    "        \n",
    "        phi_hat_d = beta*(B*alpha_array[1] - c)/(1-beta)\n",
    "        \n",
    "            \n",
    "        # solve for continuous phi \n",
    "        n_c = 501\n",
    "        pi_c_array = np.zeros(n_c)\n",
    "        pi_c_array[1:] = np.linspace(pi_d_cf, 1, num=n_c-1) \n",
    "        pi_c_array[0] = 0\n",
    "        sbar_c_array = np.zeros(n_c)\n",
    "        alpha_c_array = np.zeros(n_c)\n",
    "        \n",
    "        for j in range(n_c):\n",
    "            pi_j = pi_c_array[j]\n",
    "            sbar_c_array[j] = pi_j*s + (1-pi_j)*(1-s)  \n",
    "            alpha_c_array[j] = pi_j*alpha_H + (1-pi_j)*alpha_L\n",
    "        \n",
    "        phi_c_hat_array = beta*(B*alpha_c_array - c)/(1-beta) # continuous stationary price \n",
    "        \n",
    "        lambda_c_array = np.zeros(n_c)\n",
    "        phi_c_array_cf = np.zeros(n_c)\n",
    "        phi_c_array_cf[:3] = 0\n",
    "        w_c_array_cf = np.zeros(n_c) \n",
    "        w_c_array_cf[:3] = 0\n",
    "        ell_c_array = np.zeros(n_c)\n",
    "        \n",
    "        for j in range(2,n_c): # j = 1 = d\n",
    "            pi_j = pi_c_array[j]\n",
    "            phi_c_hat_j = phi_c_hat_array[j]\n",
    "            lambda_c_array[j] = (- ((1- beta)/(beta*eta)) * phi_hat_d * (pi_j/pi_d_cf)\n",
    "                    *((1-pi_j)*pi_d_cf/(pi_j*(1-pi_d_cf)))**k )\n",
    "            #phi_c_array_cf[j] = phi_c_hat_j + (beta*eta/(1-beta))*lambda_c_array[j]\n",
    "            phi_c_array_cf[j] = phi_c_hat_j - phi_hat_d*(pi_j/pi_d_cf)*((1-pi_j)*pi_d_cf/(pi_j*(1-pi_d_cf)))**k                                  \n",
    "        \n",
    "        for j in range(3,n_c-1):\n",
    "            ell_c_array[j] = (eta*sbar_c_array[j]*phi_c_array_cf[j+1] \n",
    "                                + eta*(1-sbar_c_array[j])*phi_c_array_cf[j-1] \n",
    "                                +(1-eta)*phi_c_array_cf[j] - c)\n",
    "        ell_c_array[-1] = ell_c_array[-2]      \n",
    "        for j in range(3,n_c-1):    \n",
    "            w_c_array_cf[j] = ((1-chi-zeta)*alpha_hat*S(ell_c_array[j])\n",
    "                                + zeta*alpha_hat*(sbar_c_array[j]*S(ell_c_array[j+1])+(1-sbar_c_array[j])*S(ell_c_array[j-1]))\n",
    "                               + chi*(sbar_c_array[j]*alpha_H*S(ell_c_array[j+1])+(1-sbar_c_array[j])*alpha_L*S(ell_c_array[j-1]))\n",
    "                                 - c)\n",
    "        w_c_array_cf[-1]=w_c_array_cf[-2]\n",
    "            \n",
    "        # closed form solution for phi              \n",
    "        lambda_array = np.zeros(n)\n",
    "        phi_array_cf = np.zeros(n)\n",
    "        \n",
    "        for j in range(2,n): # j = 1 = d\n",
    "            pi_j = pi_array[j]\n",
    "            phi_hat_j = phi_hat_array[j]\n",
    "            lambda_array[j] = (- ((1- beta)/(beta*eta)) * phi_hat_d * (pi_j/pi_d_cf)\n",
    "                    *((1-pi_j)*pi_d_cf/(pi_j*(1-pi_d_cf)))**k )\n",
    "            phi_array_cf[j] = phi_hat_j + (beta*eta/(1-beta))*lambda_array[j]\n",
    "                    \n",
    "        phi_array_cf[phi_array_cf < 0] = 0 \n",
    "        \n",
    "        # check incentive constraints  \n",
    "        if not c+tolerance >= eta*sbar_array[1]*phi_array_cf[2]:\n",
    "            print('Error: incentive constraint is wrong')\n",
    "        \n",
    "        d = 1 \n",
    "        \n",
    "        # solve for w\n",
    "        ell_c = np.zeros(n)\n",
    "        w_array = np.zeros(n)\n",
    "        sed_derivative = np.zeros(n)\n",
    "\n",
    "        for j in range(3,n-1):\n",
    "            ell_c[j] = (eta*sbar_array[j]*phi_array_cf[j+1] \n",
    "                        +eta*(1-sbar_array[j])*phi_array_cf[j-1]\n",
    "                        +(1-eta)*phi_array_cf[j] - c)\n",
    "        ell_c[-1]=ell_c[-2]\n",
    "        w_array[:3] = ((1-chi-zeta)*alpha_hat*S(0)\n",
    "                                + zeta*alpha_hat*(sbar_array[j]*S(0)+(1-sbar_array[j])*S(0))\n",
    "                               + chi*(sbar_array[j]*alpha_H*S(0) + (1-sbar_array[j])*alpha_L*S(0))\n",
    "                                 )\n",
    "        for j in range(3,n-1):\n",
    "            w_array[j] = max(w_array[0],((1-chi-zeta)*alpha_hat*S(ell_c[j])\n",
    "                                + zeta*alpha_hat*(sbar_array[j]*S(ell_c[j+1])+(1-sbar_array[j])*S(ell_c[j-1]))\n",
    "                               + chi*(sbar_array[j]*alpha_H*S(ell_c[j+1]) + (1-sbar_array[j])*alpha_L*S(ell_c[j-1]))\n",
    "                                 -c))\n",
    "            \n",
    "        w_array[-1] = w_array[-2]\n",
    "        avew = np.zeros(n) \n",
    "        for j in range(n):\n",
    "            avew[j] = pi_array[j]*w_array[-1] + (1-pi_array[j])*w_array[0]\n",
    "            \n",
    "        # solve for closed form welfare \n",
    "        # cannot have closed form solution for continuous welfare\n",
    "        \n",
    "        m1 = (1-beta*(1-eta))/(2*s*beta*eta) - (((1-beta*(1-eta))/(2*s*beta*eta))**2 - (1-s)/s)**0.5\n",
    "        m2 = (1-beta*(1-eta))/(2*s*beta*eta) + (((1-beta*(1-eta))/(2*s*beta*eta))**2 - (1-s)/s)**0.5\n",
    "        \n",
    "        if not m2>1>m1>0:\n",
    "            print('Error with m1 and m2')\n",
    "        \n",
    "        n1 = m1*s/(1-s)\n",
    "        n2 = m2*s/(1-s)\n",
    "        \n",
    "        if not n2>n1:\n",
    "            print('Error with n1 and n2') \n",
    "                    \n",
    "        Omega_cf_array = np.zeros(n) \n",
    "        Omega_d = w_array[0]/(1-beta)\n",
    "        Omega_cf_array[0:d+1] = Omega_d\n",
    "        Omega_cf_array[-1] = w_array[-1]/(1-beta)\n",
    "\n",
    "        for j in range(d+1, n):\n",
    "            sum_bw = 0\n",
    "            pi_j = pi_array[j]\n",
    "            for i in range(d, 1000):\n",
    "                if i<n:\n",
    "                    pi_i = pi_array[i]\n",
    "                    w_i = w_array[i]\n",
    "                else: \n",
    "                    pi_i = pi_array[-1]\n",
    "                    w_i = w_array[-1]\n",
    "                \n",
    "                if i <= j:\n",
    "                    b_ji = ((pi_j*(m1**(j-d))/(beta*eta*(1-s)*(m2-m1))) *(m2/(m1**(i-d-1))-m1/(m2**(i-d-1)))/pi_i)\n",
    "                    sum_bw = sum_bw + b_ji*w_i\n",
    "                else: \n",
    "                    b_ji = pi_j*m1*((m2**(j-d))-(m1**(j-d)))/(beta*eta*(1-s)*(m2-m1)*(m2**(i-d-1))*pi_i)\n",
    "                    sum_bw = sum_bw + b_ji*w_i\n",
    "            Omega_cf_array[j] = Omega_d + sum_bw\n",
    "                       \n",
    "        avew_c = np.zeros(n) \n",
    "        for j in range(n):\n",
    "            avew_c[j] = pi_c_array[j]*w_c_array_cf[-1] + (1-pi_c_array[j])*w_c_array_cf[0]\n",
    "                    \n",
    "        # iteration to solve for welfare    \n",
    "        Omega = np.zeros(n) \n",
    "        Omega[0:d+1] = Omega_d \n",
    "        Omega[-1] = w_array[-1]/(1-beta)\n",
    "        \n",
    "        Omega_test = np.copy(Omega)\n",
    "        Omega_array = np.copy(Omega)\n",
    "        diff = 10 \n",
    "        count = 1\n",
    "        while diff > tolerance and count <= 3000: \n",
    "            \n",
    "            for j in range(d+1, n-1): \n",
    "                sbar_j = sbar_array[j]\n",
    "                w_j = w_array[j]\n",
    "                Omega_test[j] = max((w_j + eta*sbar_j*beta*Omega[j+1] \n",
    "                               + eta*(1-sbar_j)*beta*Omega[j-1] \n",
    "                            + (1-eta)*beta*Omega[j]),Omega[0])\n",
    "                 \n",
    "            count += 1 \n",
    "            if count ==3000: \n",
    "                print(\"Count=3000 number of iterations reached at welfare \")\n",
    "            diff = np.max(np.abs(Omega - Omega_test))\n",
    "            Omega = np.copy(Omega_test)\n",
    "            \n",
    "        Omega_array = np.copy(Omega) # update welfare array \n",
    "        \n",
    "        aveOmega = np.zeros(n)\n",
    "        \n",
    "        # full disclosure, pi W_H + (1-pi) W_L\n",
    "        for j in range(n): \n",
    "            pi_j = pi_array[j]\n",
    "            aveOmega[j] = pi_j*Omega_array[-1] + (1-pi_j)*Omega_d\n",
    "            \n",
    "        return phi_c_hat_array,ell_c_array,pi_c_array, phi_c_array_cf,w_c_array_cf, pi_array, pi_d_cf, phi_array_cf, w_array, Omega_cf_array,Omega_array,aveOmega, avew, lambda_c_array, phi_hat_array\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "k 1.2984488866177444\n",
      "Count=3000 number of iterations reached at welfare \n",
      "k 1.1366477543541902\n",
      "k 1.1366477543541902\n",
      "k 1.0731464599173113\n"
     ]
    }
   ],
   "source": [
    "eta_grid = [0.2,0.6,0.8,0.99]\n",
    "\n",
    "for j in range(len(eta_grid)): \n",
    "    eta = eta_grid[j]\n",
    "    result = Usability(given_eta = eta).equm_allocation() \n",
    "\n",
    "    globals()[f\"phi_c_hat_array_{j}\"],globals()[f\"ell_c_c_{j}\"],globals()[f\"pi_c_array_{j}\"],globals()[f\"phi_c_array_cf_{j}\"],globals()[f\"w_c_array_cf_{j}\"],globals()[f\"pi_array_{j}\"], globals()[f\"pi_d_cf_{j}\"], globals()[f\"phi_array_cf_{j}\"],globals()[f\"w_array_{j}\"],globals()[f\"Omega_cf_array_{j}\"],globals()[f\"Omega_array_{j}\"], globals()[f\"aveOmega_{j}\"], globals()[f\"avew_{j}\"],globals()[f\"lambda_c_array_{j}\"],globals()[f\"phi_hat_array_{j}\"]= result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABFWUlEQVR4nO3dd1xV9f/A8ddhCCriBRcuxsW9BdRyFCo4s3JgaVbf+jrKsp2WP218G4plw6bYsuXMSk1NcGtLwJETBU0RJ1tk3vv5/XGBHChc5A7g/Xw8eiT3nHvu+x6v9835nM/n/daUUgghhBD2xsHWAQghhBAlkQQlhBDCLkmCEkIIYZckQQkhhLBLkqCEEELYJSdbB1CS+vXrK19fX1uHIYQQwgpiYmIuKKUaXP24XSYoX19foqOjbR2GEEIIK9A07Z+SHpchPiGEEHZJEpQQQgi7JAlKCCGEXZIEJYQQwi5JghJCCGGXJEEJIYSwS3Y5zbwscnNzSUlJITMzE4PBYOtwRDk4OjpSp04dPD09cXFxsXU4Qgg7UykTVG5uLidOnMDDwwNfX1+cnZ3RNM3WYQkzKKXIz88nIyODEydO4O3tLUlKCHGFSjnEl5KSgoeHB/Xr16dGjRqSnCohTdOoUaMG9evXx8PDg5SUFFuHJIQwQ2rKBXb9utCir1EpE1RmZibu7u62DkNUEHd3dzIzM20dhhCiFEop9iYkMSKkGyO6NaPr709w+vghi71epRziMxgMODs72zoMUUGcnZ3lPqIQdix69x5ee/sjBnVwZ1T2EkJcLvBXMy+ODPuSFj6tLfa6lTJBATKsV4XI36UQ9ic+Pp5MrSYr/z7P7m9msv7n1cyqWwtD61t4cN4SJvv3tHgMlTZBCSGEqHgFBiMLVm5l8oi+9BsyhAVB8Uxrk0Ryj1tpPOINNP9gq8UiCUoIIao5o9HI8JGjMOiak9rqTlpl7OCVoU2Z0HorDRu0xynkHZq0HgxWHu2olJMkRMWbNm0a/v7++Pv7ExERccN9ExISCA0NxcPDo0z7CyHsz59//slnn33GkbOZzPh5P5sTMkn85whfGv+PL13m8tKABjT5z5c4Tf4N2gyxenICuYISmJJTVFQUMTExAAQGBuLp6cmoUaNK3D80NJT58+cTEhJCVFQUYWFhAEycONFqMQshzPfPP//g4+OD0aiYNW8+61auoNHh+nRzOcnOcUZaZe2Amk1h0Dy0LmPB0baT0eQKShAREcGyZcvQ6XTodDqmTZvGrFmzStw3NjYWvV5PSEgIACEhIbz44ouEh4dbM2QhhJm++eYbfH19mbVoA/3f2UKMZ3/ueOoNNnl/wXKnGbRSx2HgLJgSC4EP2jw5QRW7gnp11X4OJGXYNIZ2Tdx5eVh7m8ZgjtjYWNLS0tDr9cWP6fV6YmNjS9w/ICCA+fPnX/GYXq8nISHBonEKIcxz+vRpJk2axGOPPUbboD7sVd54hYzn4z/OMdDflYVd19P8xM9oF92h3wzo8Si4uNk67CtUqQRVmUVERBATE0N4eDhRUVEAREZGMmnSJAICAiz2uikpKeh0uiseK0pWaWlp12y7fHuRnTt3WjRGIUTplFJs27aNgoIC+vXrh6enJ4eOJBD+006ObS7AUdN48MFxPO26mkYHF0I60PNx6P0M1PK0dfglqlIJqjJduVwuKiqK0aNHExMTQ1hYGJGRkcXbZs2axbJly0p83rRp00hLS7vhsXU6HS+++GKJiQYo8fmenqYPa0nJq6TnFw0RCiGsLyUlBU9PTzRN44knnkDn4UFOg7ZEbE0g765wUmo580RQQ8Y7r6dO9IeQmwFdxkLwi6Brbuvwb6hKJajKSq/Xo9PpiI6OZsGCBcWPJyQkFCeLklTEfR+dTnfdJHej1y7Sv39/FixYUHxPSghhPVOnTuXrr78mMTGRXAPcM/VtViXkMvm7WHzq1eK1Ya0Z7bwNl21PQOZpaDUY+r8EjdrZOvQykQRlB4qGzBISEq4YKisa4rOkkpJQ0f2k0q6eQkNDmTRp0nVn+wkhKtbu3bt5+eWX+fzzz6lfvz5DhgzBTVeP8F/2sWTXWTJyCgj08eClu3wJdYzBceMYuBAHzbrDqC/Ax/LVHyqSJCg7ERUVRVBQ0DWP3WjorCKG+AICAtDpdMTGxhYnx6ioqFKviEJDQwkLC5Op5UJYkNFoZOvWrXh7e6PX69E0jZiYGI4cOUJygQurz3vwc1ZHCv44xcB2Xky4zY9ADkHkg5D4F9RvBfd8B22G2mQd082SBGUnIiMjCQ0NLf65KGHc6CqmoqZ2T5w4kWnTphEZGUlCQgKzZs26YqgxLS3tiqu7sLAwSU5CWFBOTg6urq6kpqYyYMAAnnrqKebMmUOnTp34fkMsn+w4ztaft1LT2ZEx3b15uJcfvoZ/YMNkiFsLdRrDsHnQ5T5wrLxf85U38iomISHhiuG8JUuWWO2+Tnh4eHEliaKfLx+2i4iIYP78+cTHxxMbG8vy5ctZvnz5NcOPSimrxCtEVTZs2DCcnJz48ccfqVevHpGRkXTuGsiK2EQitiZw6Ewm9d1ceG5AK+7r4YNHwTnYNBX2fA816kD/l6HHI1Cjlq3fyk3T7PFLJSgoSEVHR193+8GDB2nbtq0VI7K+wMBAwsPDq83kg+rwdypESaKiovjxxx/56KOPAPjggw9wcHDgscceIzMnn0V/neCL7cc5k5FDy4ZuTOij566uTXDJS4ft78Kf8wEF3SdCn2ftdsr4jWiaFqOUCrr6cbmCslOxsbHVJjkJUZ0YDAa2bNlCr169cHFx4fDhw6xatYpXXnmFBg0aMGXKFM5n5vLWr4f4+vd/yMwp4FZ9PWaN7EhwqwZoBTnwxzxTcsrJgM5joO+LoPO29VurcJKg7FBsbKzMjBOiijEYDDg6OrJx40YGDBjAihUrGD58OBMmTODRRx/FwcGBE8mXiNgWz9LoRPINRga19+KR2/3p3FwHhgLY9Q1smgWZSdByIIS8DI0q5/rPspAEZYcCAgJk4asQVcTFixfp06cP48aN49lnn6Vv374sXbqUQYMGAVCjRg32J6Xz6ZYEftmbhJODAyMCmjLhNj3+DdxAKTi0BqJegQuHoVk3GPkZ+Pay7RuzAklQQghRwZYvX865c+eYPHkybm5uBAYG0ry5qWqDk5MTYWFhKKX4Lf4Cn2yOZ9uRC7i5ODGhj56He/vRyN3VdKDEGFg/A078BvVawj3fQps7KuWU8fKQBCWEEDfJYDCwe/duAgMDAVixYgWHDx/m0UcfRdM0Pvvss+J9jUbF+gNn+GRzPHsS06nvVoPnB7Zm3C0+1K1ZWEE85Rhs+B/sXwG1G8DQdyDgAbuoMG5NkqCEEKKclFJomsbs2bN56aWXSExMpHHjxnz88ce4u7ujXXalk1tg4Kddp5i/JYGEC1l4e9bi9bs7MCqwGa7OjqadLqXA1rfgrwWmZHTbVOj1BLjUsdE7tC1JUEIIYab9+/czduxYIiIi6NGjB2PHjqVVq1bFpcMuX2BfNFX88+3HOJuRS/sm7nwwpiuDO3jh5FjYki8/B/6aD1vnQl4mdB0HwdPBvbEN3p39kAQlhBClKCgoYNGiRTRr1oy+ffvSvHlzdDodubm5APj5+eHn53fFc1Ky8vhi+zG+/v04GTkF9PSvx9thnendov6/V1ZGI/y9DDa+BuknoeUACHm10hRztTRJUEIIUQKDwcDJkyfx9fXFwcGBGTNm0LdvX/r27Yu7uztbtmwp8XlnM3JYsDWB7/48QU6BgYHtvHg0uHCq+OUStkDkTDi9Bxp3hrs+Av3tln9jlYgkKCGEKEFYWBj79+/n0KFDODg4sH37dpo1a3bd/U+mXOLTLfEsi07EoBR3dm7C5GB/Wja66v7R2QMQ+RIcjYS6zWHEAugwChwcLPyOKh85IwKguBafv78/ERERN72/uccTwtbWrVtHz549uXTpEgCTJ0/mzTffLK4x2bx58ysmPRQ5eu4izyzdTfDbm1kWncjIwGZsejaYd+/pcmVyyjgNPz8On/aCk39B6P/g8WjoNFqS03XIFZRg2rRpREVFERMTA5jqAHp6el63msW0adNISEggPj6etLQ0/Pz80Ov1xaWZStsuhD3Iyspi2bJl9O/fn+bNm+Pq6opSiqSkJFq0aFHq53V/UjofbTrK2n1ncHFy4MFbfZl4mx6vuq5X7pibCTvmwe8fgiEfejwKtz1XKWvmWZ1Syu7+CwwMVDdy4MCBG24X5tHpdCo+Pr745/nz56uAgIAS901NTVXAFftPnTpVhYSElGn79cjfqbCG/Px8lZycrJRS6tixYwpQ7733nlnHiD6eoh768i/lM2216vDSOhW+9qA6n5lz7Y4FeUr9tUCpOf5Kveyu1NL/KJWcUBFvo8oBolUJucDqV1CapoUrpaZZ+3VFyWJjY0lLSyvu6gumDr+xsbEl7l9UZf7y/bt161Y8jFfadiFsxWg00q5dO/r06cPnn3+Or68ve/bsoWPHjqU+VynFb/HJfLjxKL8nJONRy5lnQ1vxQE/ffxfX/rszHF4DkS9D8hHw7gljlkCzQAu9s6rLqglK07QQQF/qjuW19gU487fFDl8mXh1h8GyznxYREUFMTAzh4eFERUUB/7Z8v7wNfEVLSUm5piliUXJJS0srte170f436uxb2nYhLOXLL7/kt99+Y8GCBTg4OPDUU0/h6+tbvL1Tp043fL5Sio2HzvHBxqPsPplGgzouzBjaljHdvantUsLXZ2I0rJ9pKk1UvxXcuwhaD642pYkqmtUSlKZpeiDBWq9XmURFRTF69GhiYmIICwsjMjKyeNusWbOuWzi2Ilq+l/T8osWGJSWvorb0lyevoqumtLS0UreXJeEJUV4ZGRmsXLmSsWPH4uDgwKlTpzh06BC5ubm4uLgwefLkMh3HYFSs3XeajzbFc/B0Bk11NXnt7g6EXV714XIpCYWliX68rDTRg5W6m61dKGnczxL/ASGF/192ne0TgWgg2tvb+4bjlVXtfkXR/ZqAgAAVExNT/Hh4eLiaOHGiRV87MjJSmT4G/yq6j5Samlric0JCQorvKcXHxyu9Xn/FMUrbXpKq9ncqrCcvL0/l5uYqpZT67rvvFKC2b9+ulFLKYDCYdawCg1H9tCtR9Z+7WflMW636vr1JLYs+qfIKrnOcrGSl1kxT6tV6Sr3updTGN5TKybip91MdcZ17UFaZ26hpWohSKupG+yilIpRSQUqpoAYNGlgjLLtRNKSWkJBwxXBeZGQkoaGhFn3toqulyyUkmC50r3e1U3RF5+/vX9wu/vK4S9suREU5efIkTZs25bvvvgPg7rvv5vfff6dnz54AOJRx+naBwciPuxIJfXcLTy7ejYMG88Z0JfLp2xkV2Axnx6uOk58DO96H97uYShR1GQNTYqHv9GpbN88SrHX9mVJ4/0kH6DVNC1BKlXwXvpqKiooqHh67/LEb9YWqiCG+gIAAdDodsbGxxUkkKirqhlNsdTrdFcOQYWFhTJo0qczbhSgvpRRz5szB3d2dRx99lGbNmjF69GhatWoFQK1atbjlllvKfLwCg5Gfdifx0aajHLuQRRuvOnx8XwCD2nvh4FDCfSOlTBXGo16BtBOm0kSh/4OGbSvoHYorlHRZZan/MA3jxQMBN9qvOk4znzp1qgoPDy/+OSYm5rpTvS3x2pcPyel0OrVs2bLi7ampqVcMPcbExBQP/4WHh18TZ2nbS1IV/05FxUhLS1MbN24s/nngwIFq3LhxN3XMvAKDWrLzhLptzkblM221GvzeVrX279PKYDBe/0kn/lRqQYhpyvjHPZU6uvH6+wqzYMshvsuSYYRSyl/J1dM1EhISrlgYu2TJEqstbA0PDycgIAB/f39CQ0MJDw+/IpaIiAjCwsKKf46KisLPzw8PDw/i4+OLF/iWdbsQpSkoKCj+8/Tp0xk6dCiZmZkArFy5km+++aZcx80rMLL4rxP0m7uZqcv3UsfViYj7A/nlid4M6nCdq6bU47DsP/B5qOmq6c4PYdJW8O9brhhE2Wmm5GVfgoKCVNHMr5IcPHiQtm2r9iV1YGAg4eHh1ab6QnX4OxVls2HDBu6991527NhBq1atOHLkCKmpqXTr1q3EUkNlkVdgZFnMST7eFM+ptGw6NavLk/1b0q9Nw+sfMzsNtr0Nf84HzdHUl6nnE+DiVv43J0qkaVqMUiro6sdlDqSdio2NrTbJSVRvmZmZvP/++wQHB9O7d2/atGlD3759MRgMALRs2bLcx84tMLA0OpFPNh0lKT2HLs11vD68A8GtGlw/MRnyIfoL2DwbslOhy1joNwPcm5Q7DlE+kqDsUGxs7HXr4AlRFaSmpnL69GnatWtHjRo1ePfdd1FK0bt3b5o2bcrSpUtv6vg5+QaW7DzJJ5vjOZORQ6CPB7NHdqJPy/rXT0zFFSBeguSj4HcbDHgDGt94Ma+wHElQdiggIOCGs/eEqIyUUsXJITQ0FGdnZ37//XdcXFw4duwY7u7uN/0aOfkGFv11gk82x3MuM5fuvp7MHd2Znv71bjw8mLQb1s+A49tMFSDGLIFWA6UChI1JghJCWNyCBQuYN28eu3btwsnJiTlz5lC3bt3i7TebnIqG8j7aeJQzGTn08PPkvXu7cKu+lMSUfsrUzXbPYlN18SFvQ+B/wNH5+s8RViMJSghR4U6fPk1ERASPP/449erVw8vLiw4dOpCenk69evXo169fhbxOvsHI8phEPtx4lFNp2QT5ePDOPZ3p6V//xk/MvQg73oPfPgRlhF5PQp9nwLXujZ8nrEoSlBCiQqSkpFBQUEDDhg05c+YMr776Kp06dWL48OEMGzaMYcOGVdhrFRiMrNh1ig82HuFkSjZdmuuYPbIjvVvc4B4TgNEAu76Fja9D1jnoMBL6vwwePhUWm6g4kqCEEDctKysLb29vJk+ezJw5c+jSpQuJiYk0aVKxM98MRsXKPad4P+oIx5Mv0alZXf53Vymz8ooc3WCqNH5uPzTvAWMWQbNrZjYLOyIJSghRLjNmzODMmTN89tln1K5dm/fee48ePXoAoGlahSYno1Gx+u/TvB8VR/z5LNo2dmfBA0GEtL3BOqYi5w6aJkAcjQIPXwhbCO3ukgkQlYAkKCFEmRw7doxffvmFxx9/vPgxo9FYPDtv/PjxFf6aRqNi3f4zvBcVR9zZi7RuVIdPxwUwoN11qj5c7uI52PQGxH5tKuA64A3oPgGcXCo8TmEZkqCEENeVnJxMnTp1qFGjBqtXr+bJJ59k6NCh+Pn58frrr1vsdZVSRB44y7tRRzh4OgP/BrX5YExXhnZsXHpiys+G3z+C7e9CQQ50nwS3TzXN0hOViiQoIUSJoqOj6dmzJ8uWLeOuu+7igQceYMSIETRt2tRir6mUYtPhc7wbeYS/T6XjV782793ThWGdm+BYWmIyGuHvZabGgRmJ0OYOU6Xxev4Wi1dYliQoOzVnzhwCAgKk3JGwmvz8fB577DG6devGhAkT6Ny5M88//zzt2rUDoG7dulesXapISim2H73A3PVx7D6ZRnPPmrw1qhPDuzbF6epeTCU5vgPW/x8k7YLGXWBEBPj2skiswnokQdmhiIgIQkJCiI6OLrU3kxA3Iy4ujri4OO644w6cnZ2Ji4srvkJydnbmjTfesHgMMf+k8Navh/kjIYWmuprMHtGRkSU1CSxJSoKpNNHBVeDeFIZHQMcwKGOjQmHfJEHZoZCQEPR6PQEBAcTGSmcSUbEyMzOpU8fU9XXmzJls3ryZpKQkHB0d2bRpU7krhptrf1I6c9fHsfHQOeq7ufDqne25t3tzXJwcS39yTjpsfRv+/BQcnKHvDLj1MahRy/KBC6uRBGWHilrAA9IqXVSohQsX8sgjj3D8+HEaNWrEG2+8Qa1atXB0NCUFaySn+PMXeTcyjtV7T+Pu6sTUQa35T09fatUow9eRoQBiF8KmN+FSMnS5r7DSeGOLxy2sTxKUHVu+fLlUNRc3JSkpiZdeeonJkycTEBBAjx49mDJlSlGHa1q0aGG1WBJTLzFvwxGWxyTi6uzI431bMOE2PXVrlrHuXfxG+PX/4NwB8OkFA9+EJl0sGrOwLUlQQlQxhw4dIjc3l86dO1OrVi1+/vlngoODCQgIoE2bNsyZM8eq8ZzPzOWjTUf5/s8TAPynpx+T+/pT362M65EuHDEttI1bZ1poO/obaDtMFtpWA5Kg7FRaWho6nc7WYYhKIj8/H2dnZ4xGI6GhoXTu3JnVq1ej0+lISkrC2dn61bnTL+Uzf2s8X+44Tp7BSFhgM6b0b0lTXc2yHeBSCmyZAzsXgHMt05TxHo/IQttqpMpNdQkODuarr74CTP9og4OD+fbbbwG4dOkSwcHBLFmyBID09HSCg4NZsWIFABcuXCA4OJhVq1YBcObMGYKDg1m3bh0AJ0+eJDg4mKioKAASEhIIDg5my5YtABw+fLjC3ofM3hNl9eKLLxIYGIhSCgcHBxYtWsRnn31WvN3aySkrt4CPNh2l95yNfLw5ntB2jYh65nZmj+xUtuRkyIc/PoV5XeGv+dD1fpgSa6o4LsmpWpErKDsTEREBQExMDGBKglOnTrVlSMLO7N+/nw8++IB3332XmjVr0rVrVwDy8vJwcXGhd+/eNokrJ9/A93+e4OPNR7lwMY+Qtg15JrQ17ZqUsdeTUnBkvek+U/IR0Aeb7jM1am/RuIUdU0rZ3X+BgYHqRg4cOHDD7ZVRamqqGjVqlIqPj1dKKbVs2TKllFKRkZFq6tSptgzNKqri32lF2r9/vzp37pxSSqmoqChVp04d9eeff9o4KpP8AoNa8tcJdeubUcpn2mp1z/zfVPTxFPMOcma/UgvvUupld6XmBSh1aK1SRqNF4hX2B4hWJeSCKjfEV1lNmDCBSZMmodfrSUtLK55qrtfrWb58uY2jE7agCmfanTx5kvbt2/P5558D0LdvX86cOUP37t1tGR5KKdbvP8Pg97cx9Ye9NKjjwrf/7cGiCbcQ6ONRtoNkXYDVT8OnvUxVIAbNhsl/QOtBMglCyBCfPUhLSyMqKoply5YBpvtPRdPLExISrlgXJaqH0aNHU69ePT755BOaN2/OokWL6Nu3LwAODg7UqmXbBak7j6cwe+0hYv5JRV+/Nh/fF8DgDl5lX0dVkGtaZLv1bci/BN0nwu3TpKCruIIkKDtwoyQ0f/58Jk2aZOWIhLX9/vvvbN26lWnTpgGmK+fL697de++9tgrtCofPZPLWr4eIOniOhnVceHN4R8KCyliWCEz3mQ6ugsiZkHocWg6EAa9Dg1YWjVtUTpKg7MD1qkVERESg1+tlsW4VtX//ftq1a4emafz666+8//77PProo7i7uzN79mxbh3eFU2nZvBsZxw+xibi5OPH8wNY83MuPmjXKUJaoSNJu0wSIf7ZDw3Zw/4/g389iMYvKTysa57YnQUFBKjo6+rrbDx48SNu2ba0YkeUlJCSwfPly9Ho9CQkJ6HQ69Hp9tZlqXhX/Tm/kp59+Yvjw4Wzbto3evXuTnp5OjRo1qFmzjGuErCQ1K4+PNx9l4e//gIIHe/owObgFHrVrlP0gmWdgw2uw+zvTEF7f/4OAB8FRfj8WJpqmxSilgq5+XD4hdkKv1zN16lSWL18u08qroJSUFMaMGcP999/PuHHjCAkJYd68ecVJ2VJtLMrrUl4BX+44zqeb47mYV8DIgGY8Hdqq7ItsobBx4Iew7V0w5EHPx6HPc1BTZ7G4RdUiCUoIC4mKiiIzM5Phw4fj4eFBQUEBBoMBADc3N6ZMmWLjCK9VYDCyNDqR96LiOJeZS0jbhjw/sA2tveqU/SBKwb4fIOoVSD8pjQNFuUmCEqICJSYm0qxZMwBmz55Namoqw4cPR9M0NmzYYOPork8pxbp9Z3jr18MkXMgi0MeDj+4LoJuvmbPqEmNg3QuQ+Bd4dYK7PwG/PpYJWlR5kqDsjEyIqLxeeeUV3nrrLc6ePYubmxtffPEFDRs2tHVYpfrrWApvrjnI7pNptGzoRsT9gYS2a2Re642M06Yrpr2Lwa0R3PkhdBkLDmZMohDiKpKghCinffv28fTTT/PJJ5/QokUL7rrrLho0aFD8xe7t7W3jCG8s4fxFZq89xPoDZ2nk7sKckZ0YEVDGFutF8nMK7zO9A8Z86P009HkWXMwYEhTiOiptglJKWa3zp7Ase5xJWhKDwUBkZCSNGzemc+fO6HQ6jh8/zsmTJ2nRogVdu3Ytrotnz5Iv5vL+hiN8/+cJXJwceG5AK/7bW2/elHGl4OBKUxuMtBOm+0wDXgdPP8sFLqqdSpmgatSoQXZ2ts1X04uKkZ2djYuLfVapVkqRkZFB3bp1yc/P59577yUsLIwFCxbQrFkz4uLiKs0vStl5Br7YcYxPNseTnW9gTPfmPNm/FQ3qmHnuT++FdS8WrmdqDw+uAr/bLBO0qNYqZYKqX78+iYmJ1K9fnzp16uDk5FRpviSEiVKKgoICMjMzuXDhAo0aNbJ1SCUKCwvjzJkzbN++HVdXVzZt2kS7du2Kt1eGz53RqFix6xRz1x/mdHoOoe0aMW1QG1o0dDPvQFkXYONrELMQanrA0HdkPZOwqEr5yapbty4uLi6cP3+e5ORkCgoKbB2SKAcnJydcXV3x9vbG1dXV1uEAsGHDBiIiIvj+++9xdHRkxIgRZGRkFA8pV4YhvMttP3KBN9cc5MDpDDo3q8t793Shh76eeQcpyIO/IkzNA/Oz4JZH4fappiQlhAVVygQF4OrqSvPmzW0dhqjkCgoKWL9+PT179kSn05GcnExMTAyJiYn4+PgwduxYW4dYLofPZPLmmoNsiTtPM4+avH9vF4Z1aoKDg5lXfHHr4dcXIfkotAg19WeSunnCSiplqSMhboZSqri5X0xMDEFBQSxYsIDx48djMBhwcHCoFEN3JTmbkcM76+NYFnMSNxcnpvRryQM9fXBxMnO69/nD8Ot0OBoF9VqaElOrAZYJWlR7UupICCA/P5+goCAGDx7M7NmzCQgIYO3atfTrZypa6uhYOdftZOUWMH9rAgu2JlBgNPJQLz+m9GuBrpYZNfMAslNhczjsXADOtU2JqdsEcDLzOEJUAKslKE3TiqqehiqlplnrdYVYunQphw8fZubMmTg7OzNkyBC6dOkCmCY5DBo0yLYB3gSjUbE8NpG3fj3M+cxc7ujUmKkD2+Bdz8wZroYCiP0KNr4BOWmmyQ/9ZkDt+pYIW4gysUqC0jQtAAhQSs3RNG2apml6pVSCNV5bVD/5+fns2LGD4OBgALZv387mzZuZPn06jo6OzJo1y7YBVpA/E5J57ZcD7DuVQVdvHfPvDyTAuxwTFxK2mKaNn9sPvn1g0Czw6ljxAQthJqskKKVULBCraZoOSJDkJCqaUgqlFA4ODsyfP58pU6YU91sKDw/H1dW10t5XutqJ5EvMWnuQtfvO0KSuK+/f24U7Ozcx//2lHDMttD20GnTeMPobaDtMWq0Lu2Hte1BBQHxJGzRNmwhMBPsvESPsS3x8PEOHDuWtt95i2LBh3HPPPfj4+NCiRQsAu+uxVF6ZOfl8uOkoX24/jqODxrOhrRjfx8wKEAC5mbBtLvz+ETg4Q/+X4JbHwNk+pvoLUcSqCUopFaVpWpimaaOUUsuv2hYBRIBpFp814xKVi9Fo5PPPP6devXqMGDECb29vWrVqVVxZpEGDBgwbNszGUVYcg1GxZOdJ5q4/THJWHiMDmjF1UGsauZuZUIxG2LMINrwKF89C5zHQ/2Vwb2yZwIW4Sda6BxUOxBcmoTTAzBr+orrLy8sjLi6ODh064ODgUFygdcSIETg7O7Ny5Upbh2gRO45e4LXVBzh0JpNuvh58+VA3OjXTmX+gE3/CummQtAuadYN7F0GzwAqPV4iKZK0rqPmAvnAmn64wUQlRZhMnTmTNmjWcOnUKZ2dnIiMj8fSsur/nJJy/yJtrDhF18CzNPGry8X0BDO7gZf59pvREiHwZ9i2HOk1gxALoGCb3mUSlYK1JEglA0cSIKGu8pqjcNmzYwBNPPMHmzZtp0KABjz32GGFhYTg4mFpB1KtnZrmeSiL9Uj7zNh5h4W/HcXFyYOqg1jzcyw9XZzPvM+Vnw28fmNpgoOC2qdD7KahR2xJhC2ERslBX2IWsrCy+/fZb+vTpQ7t27WjUqBENGzbk/PnzNGjQgG7dutk6RIsyGBXf/3WCd9YfJi07n3uCmvPMgFY0rGPmfSal4OAqWP9/pjYY7e6GAa+ZZukJUclIghI2k5eXR3JyMo0bNyYvL48nnniCmTNn0q5dOzp06MCmTZtsHaJV/JmQzCurDnDwdAa36D2ZeUc72jepa/6Bzh2EtdPg2JbCNhirpd26qNQkQQmb6datGz4+PqxcuRIPDw8OHjyIn1/1aXiXlJbNrLWHWLUniaa6m7jPlJ0Km2fDXwtMnWyHvA2BD0kbDFHpySdYWM3ChQv54Ycf+Pnnn9E0jRdeeOGKiQ56vd6G0VlPTr6BBVsT+HhzPEaleLJ/Sx653d/89UxGA+z6Bjb8z5SkAh8ylSeqVXUnj4jqRRKUsJjU1FSWLl3Kgw8+iKurK7m5uWRlZZGZmYm7uztjxoyxdYhWpZRi/YGzvP7LAU6mZDO4gxfTh7SluWc5OkOf+APWToXTe8C7JwwOh8adKj5oIWxI2m2ICpWbm0t+fj5ubm5ERkYyYMAAfvnlF4YMGWLr0Gzq6LlMXl11gG1HLtCqkRsvD2tPrxblKMSakWSaNv73UnBvapoA0X6ETBsXlZq02xAWl5KSQosWLXjhhReYOnUq/fr1Y9euXXTu3NnWodlMRk4+70eZpo3XrOHIy8PaMe4WH5wdHcw7UEEu/P4hbJ0LxgK47Xno/bRMGxdVmtkJStM0XyBEKfVZ4c9dMBWAzajY0ERlEB4eTk5ODi+//DKenp488cQT9OrVCzD1Vipqa1HdGI2K5TGJzPn1EMlZedzbrTnPDWhNPTcX8w6kFMStM1UbTz0Gbe6AAa+DZ/WZTCKqL7MSlKZp44FHgLrAZ4UP+wMvAvdUbGjCHl24cIFt27YxfPhwAA4cOEBOTk7x9ldeecVGkdmPPSfTeOnnfexJTCfAW8eX/+lOx2blmDZ+4Qise8HU1bZ+a7j/R/DvV/EBC2GnzLoHpWnaEUwVyWOUUi0uezxZKVVhS/vlHpR9ycnJwcXFBU3TePnll3nttddISkrCy8sLo9FYXN2hukvNymPOr4dZvPMEDdxcmD6kLXd1KUcbjJwM2DoH/vgEnGtB8IvQfQI4OlsmcCFsrKLuQXkqpdJL+Acnd2irqN9++40hQ4awZs0aevbsyaRJkxg5ciReXl4AkpwwDectjT5J+LpDZOQU8HAvP54KaUkdVzMTSlG18ahXIOs8dB1nqjbu1sAicQth78xNUBs0TRsBFF92aZq2BKmvV2VkZ2fzxhtv0L17d+688046duzIiBEjqFvXNETVpEkTmjRpYuMo7ce+U+nM/Hkfu06k0d3Xk//d3Z42Xu7mHygxBtY+D6diTNXGxy6BpgEVH7AQlYi5CWoCsAHw1zTtV0zDfQlA/4oOTFjP+fPnOXbsGN27d8fFxYXFixejlOLOO++kTp06fPHFF7YO0e6kZ+fzzvrDfPPHP3jWrsHcsM6MCGhq/nBe5llTf6bd34FbIxg+HzqOBrkyFcK8BKWUSgeCNE3rD+iBOUqpDRaJTFiUwWDA0dFUuWDs2LEcP36cuLg4HBwc2LdvH66u0l21JEopVsSeYtbag6Rk5XH/LT48M6A1dWuaOZxXkAd/zYfN4VCQA72eNE0dd6ljmcCFqITMncXnC1CYlDYUPtYP0zTz4xUdnLCMr7/+munTp3P48GFq167Nm2++Sc2aNYt/+5fkVLJDZzKY+dM+dh5PpUtzHV891J0OTcsxO+9oFKx9AZKPQMuBMPBNqN+i9OcJUc2YO8Q3HwgHjl/2mEfhYzLN3E6dOXOGefPmMWHCBPz8/GjRogUhISFkZGRQu3btKt/K4mZl5uTzbuQRFv5+HHdXJ8JHdiQssDkODmYO56Ueh3XT4fAv4OkPY5dCq4EWiVmIqsDcBBWklNp4+QNKqR80TZMOuXbm7Nmz5OTk4OPjQ15eHm+//Tbt27fHz8+Pnj170rNnT1uHaPeUUqzee5rXVh/g/MVc7u3mzdSBrfGoXcO8A+Vnw/b3YPu74OAEIa/ALZPBycxFu0JUM+YmKE3TtDpKqcyrH6+ogET5KaXQNI2CggLat2/PHXfcwVdffYW3tzdnz57Fw8PD1iFWGidTLjHjp31siTtPh6buRDwQRJfmOvMOohQcXmNabJt2AjqMhNDXoG5Ti8QsRFVjboJahmk4b3LRA5qmfQIsrcighPleeeUVoqOjWb16NU5OTsyfP5+2bdsWb5fkVDb5BiOfbTvG+xvicNQ0XrqjHQ/c6oOTubXzLhyFddNM95satJXmgUKUg7mz+CZpmhajaVoypunlemSauU0cPXqU77//nhkzZuDg4ICnpydeXl7Fs/NGjhxp6xArnZh/Uvm/H//m0JlMBrRrxCt3tqeJrqZ5B8nLgq1vwW8fgnNNGDhLqkAIUU7larehaVoI4Idp9l6FTzOXUkclO336NHXq1MHNzY3Fixczbtw4YmJiqnW18IqQnp3PnHWH+P6vE3i5u/LKne0Z2N7LvIMoBft/hPUzIOMUdB5rutdUp5FFYhaiKqnQdhtKKakcYWWHDx+mXbt2zJ8/n/Hjx3P33XeTlJREw4YNbR1apVU0CeLVVQdIycrloZ5+PDOgFW4uZv6zOHfIVAXi2Fbw6gijvgDvWywTtBDVyA3/JWqaNguYX7TGqbCaeYmK2m+IiqGUYvz48fj7+zN9+nRatWpFeHg4ffv2BUxrlWS9UvldPgmiY9O6fPVQN/PXNOVkwJZw+PNTqOEGQ+ea2q47mNm6XQhRotJ+VQwDIvl33dMj19lP8W/7DVFO+/fvZ9euXYwbNw5N07h06RLZ2dkAaJrGc889Z+MIK7+rJ0G8PKwdD9zqi6M5a5qUgr1LIXImXDwHAQ+YirrWrrCC/kIIzG+34W6NxoTV6R7U+fPnadDAVK368ccf5+uvv+bcuXNydWQBe06mMe2HvRw6k8nA9qZJEI3rmjkJ4vReWDsVTvwOTQNhyFum/wshyu1696DMrUh5TNM0KRZWQZYtW4aXlxeHDh0CYPr06cTHx0tyqmDZeQbe+OUAwz/eQeqlPCLuD2T+/UHmJafsVPjlOYi4HS7EwZ0fwH+jJDkJYUHmTpJYDswBHrVALFXe+fPneeaZZ3jggQcIDQ2lT58+TJ8+/YpWFqJi/Xb0Ai+s+JsTKZcY28ObFwa3wd2cPk1GI+z+1tSjKTsVuo2HvtOhpqwrE8LSzE1Q64EFmqbpMd2bSivaIJMkSrZ7924uXrxI7969qVu3Ln/88QfBwcEAeHl58dprr9k2wCoqPTufN385yJLok/jWq8Xiibdwi97Me0SnYkxXTUmx4H2raTjPq6NlAhZCXMPce1DXuzGklFIVVnG0st+DysrKonbt2gB06dKF2rVrs2PHDgBpkW4F6/ad4aWf95Gclcf4Pn48HdIKV2czZtZlJZt6NMV+DW4NTeWJOo0Gc3s9CSHKpELWQZV0AHGlV199lYiICP755x+cnJxYuHAhzZo1K94uyclyzmXm8MrK/az5+wxtG7vz+YPd6NjMjKnjRgPEfAkbXoPcTLj1Mbh9GriWo0OuEOKmlWuhbmEPKD0Qr5TaVLEhVS4HDhzgzTff5L333qN+/fr06tULpRS5ubk4OTlJlQcrUEqxPCaR1385SHa+gecHtmbibXqczamflxgDvzwDp3eDbx/TcF7DtqU+TQhhOeY2LOwKRGHqAZUA6DVNOwoMqC4NC5VS7Nq1C09PT3x9fSkoKGDNmjWMHz+e4OBgQkJCCAkJsXWY1captGxe+GEv245cIMjHg9kjO9GioVvZD3ApxTScF7MQ6niZqkC0HyHDeULYAXPvQR0FYpVSoy97bBlQRyk1qKKCssd7UAUFBTg5OZGamkqjRo2YMmUKc+fOBSAvL48aNczsESRuilKKJTtP8vovBzEqxQuD2zCuh0/ZmwgajbDrG9PsvJx0uOVRCH5BWq4LYQMVVYvP8/LkVGgCkFLuyCqB++67j+zsbFasWIGHhwcrV66ke/fuxdslOVnX6fRspv3wN1vjznOrvh5zRnWiuWctMw6wB355FhJ3mmbnDZ0LjdpbLmAhRLmYm6CiSmhYqDAN+1UZv/32Gz/++CNvvfUWAIGBgeTm5hZvHzSowi4WhRmUUiyLSeS1VQcoMCr+d1d7866actJh4xuwcwHUqgd3fwqd75XhPCHslLkJKgXYqGnaksseCwVSNU0rLhSnlHq7IoKzFqUUsbGxdOjQARcXF3bt2sUXX3zBs88+i5eXF88884ytQ6z2zqTn8OKKvWw6fJ7uvp68FdYJn3q1y/bkotp562fApQsQ9F/oNwNq6iwasxDi5lTUOqjL3fSaKGvdgypqkb5p0yb69evHDz/8wIgRI8jOzsbR0VGG7uyAUooVsad4ddV+8gxGpg1qw4O3+pb9quncQdNi23+2m8oSDX0HmnSxaMxCCPPIOqjLZGdnM3jwYO68806eeeYZ+vTpw+eff17cyqJmTTMLiAqLOJeRw/Qf/ybq4DmCfDx4K6wzfvXLeNWUmwmbZ5taYbjUgWHvQ9cHQNahCVFplGsdVGW0bt06kpKSePjhh6lZsyZNmjRBp9MB4OTkxMMPP2zbAMUVVu5JYuZP+8jJNzBjaFse6uVXtpYYRZ1tf/0/yEwqbIXxirTCEKISKlfLd7NfRNN0mBb26oFuSqlpN9q/Iob4lFLExcXRunVrAMaMGcOuXbs4ePAgmtwUt1vpl/KZ+fM+Vu5Joqu3jrfDOuPfoIzrmi4cgTXPQ8Im8OpkGs5rXmEVuIQQFlKhLd/LYTSAUipC07RumqZNVEpFWPIF586dy7Rp00hMTKRx48bMmzcPnU4nycmO7Th6geeW7eF8Zi7Phrbi0WB/nMpSDSLvEmybCzveB+daMPgt6PZf6WwrRCVnlQR1VTIqqoRuUSNGjMDDw4M6dUwLL4uaAgr7k5NvYM66w3yx4xj6BrVZMbknnZrpyvbkQ2tg7TRIPwGd7oUBr5kKvAohKj2r3oMqbNORopS6Zt2UpmkTgYkA3t7eN/1aer0evV5/08cRlrU/KZ2nFu/myLmLPHCrDy8ObkvNGmW48kk5ButegLh10KAt/GcN+PayfMBCCKux9iSJUUqpSSVtKLzKigDTPSirRiWszmBURGxN4J3Iw3jUqsFXD3UjuHUZrnwKcuG3ebD1bXBwggGvQ49HwNGMJoRCiErBaglK07RRSqk5hX8OUErFWuu1hX05mXKJZ5fu4a/jKQzu4MWbwzviUbsMa86Ob4fVT5tarre7CwbOgrpNLR+wEMImrJKgNE0LAcI1TXux8KEbzuITVddPu04x46d9aMA7ozszvGvT0ieuZF2A9TNhz/eg84Gxy6DVAKvEK4SwHWtNkogC/K3xWsI+Xcwt4KWf9rFi1ym6+Xrw7j1daOZRSoHXoorjkS9BXhb0eRb6PAc1zCgMK4SotKrNQl1hO3sT03hi0S5OpFziyf4tmdKvRenTx8/uNw3nnfwTfHqZ1jQ1bGOdgIUQdkESlLAYo1Hx+fZjzPn1EA3cXFg88Va6+3ne+El5WbAlHH7/CFzc4e5PoPMYqTguRDUkCUpYxPnMXJ5dtoetcecZ2L4R4SM7oatVykSIw2tNlSDST0LX+yH0f1CrlIQmhKiyJEGJCrcl7jzPLt1NZk4Br9/dgft6eN94IkR6ommx7aHVpjVND60Dn1utF7AQwi5JghIVJq/AyNvrDxOxNYFWjdz4bvwttPa6QQt1Q76p2vimWaCMEPIq3PqYrGkSQgCSoEQFOZWWzWPfxbL7ZBrjbvFmxtB2uDrfoCLEyZ2w+ik4uw9aDoQhb4GHj9XiFULYP0lQ4qZtPHSWZ5buwWBQfHxfAEM6Nr7+ztmpEPUqxHwF7k3gnm+hzR0yCUIIcQ1JUKLcCgxG5kbG8cnmeNo1dufj+wLwvV5DweK26/8Hl1JMQ3nBL5iaCQohRAkkQYlyOZuRw5RFu/jrWApjunvz8rAbDOklx5vWNB3bYmq7Pm4FNO5k3YCFEJWOJChhth1HL/Dk4l1k5Rp4957ODO/arOQdDfmmwq5b5oBjDRjyNgQ9LH2ahBBlIglKlJnRqPhg41He2xBHiwZuLJoQQMtG1xmiO7kTVj0J5/ZD2zth8Bxwv8G9KSGEuIokKFEm6ZfyeXLJLjYfPs/wrk15Y3gHatUo4eOTkwEb/gc7P4M6jeHe76HNUOsHLISo9CRBiVIdOpPBpG9iSErL5rW7OzDuegtvD642VYLIPA3dJ0K/GeDqbv2AhRBVgiQocUOr9ybx/LK91HF1YvHEWwj0KaH0UEaSKTEdWg2NOsA930CzIOsHK4SoUiRBiRIVGIzM+dVUFSLQx4NP7gugobvrlTsZDRD9hWldkzFfKkEIISqUJChxjZSsPKYsimXH0WTuv8WHmXe0o4bTVe0xzu6HlU/AqWjQ94U73gFPvW0CFkJUSZKgxBX2nUpn0jcxnL+Yy1ujOhEW1PzKHfKzTdPGf5sHrnVheAR0Gi2VIIQQFU4SlCi29u/TPL10N561arD8kVvp1Ex35Q4Jm2HVU5B6DLrcB6GvQe16NohUCFEdSIISKKX4cONR5kbGEeCtY/79QTSo4/LvDpdS4Nf/gz3fm4bxHlgJ+tttF7AQolqQBFXN5eQbmLp8Lyv3JDGia1PeHNHx35JFSsGBn0wz9LJToc+zcNvz4FzTpjELIaoHSVDV2LmMHCZ8Hc3eU+lMG9SGR27X/7u+KfMM/PKsaep4485SP08IYXWSoKqpfafSGb8wmoycfOaPC2RAey/TBqVg17emIT1Drqnt+i2PgaN8VIQQ1iXfOtVQ0WSIerVd+OHRnrRtXFjtIeWYqX7esS3g0wvu/ADq+ds2WCFEtSUJqpr5bFsCb6w5SIC3B/PvD6S+m4tpwe2f82Hja6A5wtB3IPAhcHAo/YBCCGEhkqCqCYNR8fovB/hyx3GGdPTindFdTJMhzh2Enx83LbhtOdC04LbuddpnCCGEFUmCqgZy8g08uXgXv+4/y397+/F/Q9riYMyHzW/D1rdMXW1HfAYdR8mCWyGE3ZAEVcWlZOUxfuFOdp1M46U72vFwbz84FWO6ajp3ADqMgsHhULu+rUMVQogrSIKqwv5JzuI/X+4kKS2bT+4LYFCruqbZeX98DG5eMGYxtB5s6zCFEKJEkqCqqL2JaTz05U6MSvH9hB4EanHw6WRIiTdNgAh91VRLTwgh7JQkqCrot/gLTFgYjUftGnzzQEf89s6F3z8CXXMpUySEqDQkQVUx6/ef4fFFu/CtV4tFgxyot3wQJB+FoP+arppc6tg6RCGEKBNJUFXIDzGJTP1hLwFNXPnGdy2uSz6Fus3hgZ9BH2zr8IQQwiySoKqIL7Yf43+rD/Bg87O8bPgQh+h4CHrYVKpIrpqEEJWQJKhKTinFe1FH+HTDfiIarSH0wnI096Zy1SSEqPQkQVViSilmrz3EX9t+ZZv75zRMP2GaoTfgNblqEkJUepKgKimlFLNW7cbzr7n84PILWs0mcM9P4N/X1qEJIUSFkARVCSml+HTJz4w4MIM2TidRAQ+iDXgdXN1tHZoQQlQYSVCVjLGggKjPp/PfpM/IddGhwpahtRpg67CEEKLCSYKqRIzJxzjx+f0MuPQ3B+v1p81/F6DVrmfrsIQQwiKk4U9loBQq+ivyP7oVz6x4Vrd4lTaPL5fkJISo0qyWoDRNG6VpWqS1Xq/KuHgOtegetNVPEp2v55uuixl635No0kxQCFHFWW2ITym1XNO0SdZ6vSrh4CpY9SSGnIu8kX8/dJ/ES3d2QJOeTUKIasBufg3XNG2ipmnRmqZFnz9/3tbh2FZOBvz4KCwZxzmHhgzMfp2srhOYOUySkxCi+rCbBKWUilBKBSmlgho0aGDrcGzn5E74tDfsXcwe/UR6XniR9p26MWtEJxwcJDkJIaoPu0lQ1Z7RYGq//sVAQLG19zfcdSCYvu2aMnd0ZxwlOQkhqhmZZm4P0hNhxUT4Zwd0GMnvbWfy3+8P0t3Pgw/GdMXZUX6PEEJUP1ZLUJqmhQBBmqaNUkott9br2r0DK2HlFDAWwN2fsq/+YMZH/IF/AzcWPBCEq7OjrSMUQgibsOYsvijAw1qvZ/fysmDdixC7EJoEwMjPOEFj/vPJb+hq1eCrh7pTt6azraMUQgibkSE+Wzi9F374L1w4Ar2fhuDpJOcoHvjkNwqMRhY/3AOvuq62jlIIIWxKEpQ1KQUxX8LaF6CWZ2HPptvJyTcw4es/OJ2ew/cTetCiobTKEEIISVDWkpsJq56CfcvBvz+MiIDa9VFKMe2HvcSeSOPj+wII9PG0daRCCGEXJEFZw5l9sOxBSEmAfjOg97NQWKroo01H+Xl3Es8NaMWQjo1tHKgQQtgPSVCWpBTEfg1rp4JrXXhgJfj1Kd685u/TvL0+jru7NOGxvi1sGKgQQtgfSVCWkpcFq5+BvYvB73YY+Rm4NSze/HdiOs8s3U2At47ZIztJCSMhhLiKJChLSEmAxePg3AEIng63PQcO/65nSsnK45FvY/CsVYP598taJyGEKIkkqIp2JAp+eBjQYNwP0KL/FZsNRsWTi3dxPjOX5Y/eSoM6LraJUwgh7JwkqIqiFGybCxtfh0bt4Z5vwdPvmt3ei4pj25ELzB7RkU7NdNaPUwghKglJUBUhNxN+etTUv6nDKLhzHtSofc1uUQfO8sHGo9wT1Jx7u3vbIFAhhKg8JEHdrOR4WDzWVBVi4Jtwy2QoYcLDP8lZPL10Nx2b1uXVu9rbIFAhhKhcJEHdjGPbYMk40Bzg/h9Bf3uJu+UVGJmyaBcOmsbH9wXIpAghhCgDSVDltetbWPUk1GsBY5eAh+91d527/jB7E9P5dFwAzT1rWS9GIYSoxCRBmctohA2vwo73QN8XRi80LcK9jm1HzjN/awJje3gzqINUihBCiLKSBGWOvEvw40TTZIjAh2DIW+B4/ZYYyRdzeWbpHlo2dGPm0HZWDFQIISo/SVBldSkFvr8HEnfecDJEEVMR2L9Jz87n64e7U7OG3HcSQghzSIIqi4wk+GYEpMSbhvTa3VXqU37cdYqog2eZMbQtbRu7WyFIIYSoWiRBlebCUfhmOGSnwH3LrztT73JnM3J4ZeV+gnw8eKjXtYt1hRBClE4S1I0k7YZvR5r+/J/V0KRrqU9RSjF9xd/kFhiZM6oTjg5SBFYIIcrDwdYB2K2kXfD1neBcCx7+tUzJCWBF7Ck2HDrH8wNbo2/gZuEghRCi6pIrqJIk7YKv7zJNH39wNXj4lOlpFy7m8uoqGdoTQoiKIFdQVytncgJ4c81BsvMNzB7ZUYb2hBDiJkmCuty5g/D13eVKTn8kJLMi9hQT+uhp0bCO5WIUQohqQhJUkfRE04QIJ1ezk1O+wcjMn/bRVFeTKf1aWjBIIYSoPuQeFJgW4X4zwtQ246G1ZiUngC+2H+PIuYsseCBIFuQKIUQFkQRVkGdql5F6DMatAK8OZj39dHo270UdIaRtI0LbNbJQkEIIUf1Iglo3DU78DiM/B78+Zj/97V/jMBgVLw+TWntCCFGRqvc9qOgvIfoL6PUkdBxl9tP3J6WzYlciD/XylTYaQghRwapvgjq7H9ZOA//+0P9ls5+ulOLNNQepW9OZyX1bWCBAIYSo3qpngsrPgR8mgKs7DJ8PDuZPbNgcd54dR5N5sn9L6ta8fssNIYQQ5VM970Ft+B+c228q/urWwOynG42K8LWH8K1Xi/t6mDfjTwghRNlUvyuoUzHwx8fQbTy0DC3XIdbuO8OhM5k8HdqKGk7V7xQKIYQ1VK9vV6MBVj8Dbo3Kdd8JTFdP72+Io0VDN+7o1KSCAxRCCFGkeiWoPYvh9G4Y+Ibp/lM5rNl3mrizF3mif0uptyeEEBZUfRJUQR5smQ2Nu0CHkeU6hMGoeD/qCC0aujG0Y+OKjU8IIcQVqk+C2v0tpJ2AfjNBK9+Vz9p9pzly7iJPytWTEEJYXPVIUErB7x9D00Bo0b+ch1BEbE1AX782Q+TqSQghLK56JKjj2yH5iGnmXjmvnnYeT2VvYjoP9/aTqychhLCC6pGgor8w9XhqP7zch1iwLQGPWs6MDGhWgYEJIYS4HqslKE3TRmmaFqJp2kRrvSYA2WlwcBV0HgPONct1iGMXsog6eJZxt/hIOw0hhLASqyQoTdNGASilogp/DrHG6wJwNAqM+eWeuQemfk/ODg7cf6tUjRBCCGuxVqmjbsCSwj8nAAFAlKVe7MyJI3h+3gMAJ4wkU5den55DsaZcx8s3KMICm9GwjmtFhimEEOIGrJWgdFf9XO/qHQqH/oqG/y5qmna4Al63PnAB0oFhN3Wgtwv/qwIKz4m4jJyTa8k5uZack2tV1DkpcXjKWgkqDfC80Q5KqQggoiJfVNO0aKVUUEUes7KTc3ItOSfXknNyLTkn17L0ObHWJImd/HsVpQcirfS6QgghKimrJCil1HJAXzg5Qlc0WUIIIYS4Hqv1g1JKzSn8ozWTU4UOGVYRck6uJefkWnJOriXn5FoWPSeaUsqSxxdCCCHKpXpUkhBCCFHpVImW74ULgdMAfeFsQLO2V0U3es+apukwTVbRA92UUtOsHqCNlPWzoGlaeHU5L2X49xOA6bNSdD+5ypPvlGsVvudJSqkSW5Fb4pxU+iuo0qpU2LSKhY2U4T2PBoKKvmysXn7KRsr6WSh8XG/F0GymjOfkxcLPiqemaVX+vJThOyUESCjcnlCYwKu8G/1yYqnv2UqfoDBVqUgo/HNRlQpztldFN3zPSqmIy37D0V+2b1VX6meh8Au4upwPKOWcFP7yslPTNH3h56Y6nJvSPifRwLKiK0ulVKw1g7NTFvmerQoJSnfVz1dXqShte1Wku+rnEt9z4ZdxSjWa9q+76ueSzou+mnwJF9Fd9fPV58S/8LEUTdPmFw4PV3W6q36+4pwopdKA+cAyINA6Idk93VU/V8j3bFVIUGncuEpFadurojTK9p5HKaUmWTgWe5LGDc6Lpmkh1ShZF0mj9M9KfOGXcgz/liOrytIo5XMCRCml/IG0ouGtai4NC3zPVoUEVVqViupYxaLU96xp2qiitWnVZQyd0s9LSmFLmFGYFpZXh/NSln8/RXSYvoiqutLOScBlw3qzqH6/AJfEIt+zlT5BXa9KhaZpkTfaXpWVdk4KHw/XNC1G07QYqsk/sDJ8VmILH/Pk2iGLKqmM/350RTe9q8OMtdLOCRChadrEwu2jq8M5geLvjaDLrxgt/T0rC3WFEELYpUp/BSWEEKJqkgQlhBDCLkmCEkIIYZckQQkhhLBLkqCEEELYJUlQQggh7JIkKCGEEHZJEpQQQgi7JAlKCCGEXZIEJYQQwi5JghLChjRNG6VpWrymaalX/RevaVq4reMTwpaqRMt3ISqjou60Sil/TdMmKqUiiiqoSxM8IeQKSghbSrmsjbau8P8hVK+OvkJclyQoIWyksAlg8ZVUIf+ix4Wo7iRBCWF704Ci/jn6G+0oRHUiCUoI2wu5/J5TUXNAIao7SVBC2FBhd9Lllz0USzXp5itEaaSjrhBCCLskV1BCCCHskiQoIYQQdkkSlBBCCLskCUoIIYRdkgQlhBDCLkmCEkIIYZckQQkhhLBLkqCEEELYpf8HtcOVKKr41y4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(pi_c_array_0,phi_c_array_cf_0,label=r'$\\eta$ = ' + str(eta_grid[0]))\n",
    "plt.plot(pi_c_array_3,phi_c_array_cf_3,label=r'$\\eta$ = ' + str(eta_grid[3]))\n",
    "\n",
    "plt.plot(pi_array_1,phi_hat_array_1,label=r'$\\hat \\phi$',linestyle=':',color='k')\n",
    "\n",
    "#plt.axvline(x=pi_d_cf_0)\n",
    "#plt.axvline(x=pi_d_cf_1)\n",
    "plt.xlabel(r'$\\pi$', fontsize=16)\n",
    "plt.ylabel('price', fontsize=16)\n",
    "plt.ylim([0,phi_c_array_cf_0[-1]*1.05])\n",
    "plt.legend(fontsize=16)\n",
    "plt.tight_layout()\n",
    "\n",
    "po = next(i for i, x in enumerate(phi_c_hat_array_0) if x > 0)\n",
    "\n",
    "# DataOut = np.column_stack((pi_c_array_0,phi_c_array_cf_0))\n",
    "# np.savetxt('usability_phi_loweta_lowalpha.dat', DataOut)\n",
    "# DataOut = np.column_stack((pi_c_array_3,phi_c_array_cf_3))\n",
    "# np.savetxt('usability_phi_higheta_lowalpha.dat', DataOut)\n",
    "# DataOut = np.column_stack((pi_c_array_0[po-1:],phi_c_hat_array_0[po-1:]))\n",
    "# np.savetxt('usability_phi_hat_lowalpha.dat', DataOut)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "8\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA/h0lEQVR4nO3deVxU973/8deXfXEZEERFAYe4L4mARhPNJmQ3zQKapU2XROiWJreLdMltb5umFtv7a29vk1ux2dpmMZCt2QVNNDGbgEncTRwENxQFBNmX7++PMyCg7DNzZobP8/HgEWY758ME5813Od+v0lojhBBCuBsfswsQQgghzkcCSgghhFuSgBJCCOGWJKCEEEK4JQkoIYQQbsnP7AIGIyIiQsfFxZldhhBCCAcoLCw8qbWO7H6/RwZUXFwcBQUFZpchhBDCAZRSJee7X7r4hBBCuCUJKCGEEG5JAkoIIYRbkoASQgjhliSghBBCuCUJKCGEEG7JI6eZCyF619jYSEVFBTU1NbS2tppdjhiGfH19GTlyJOHh4QQGBg7qGBJQQniZxsZGSktLCQsLIy4uDn9/f5RSZpclhhGtNc3NzVRXV1NaWkpMTMygQkq6+IagubWN+5/bzqeHqswuRYgOFRUVhIWFERERQUBAgISTcDmlFAEBAURERBAWFkZFRcWgjiMBNQTbiit45dOjPPLOl2aXIkSHmpoaRo0aZXYZQgAwatQoampqBvVaCagh2Lj3BADv7D3ByTONJlcjhKG1tRV/f3+zyxACAH9//0GPg0pADcGmvSewRobS0qZ5efsRs8sRooN06wl3MZTfRQmoQbKVn6H4ZC3fvCSOCydZyCk4jNba7LKEEMJrSEAN0iZ7996V08eSljiRfcdr2Hmk2uSqhBDCe0hADdLGPSeYPm4kE8NCWHbhBAL9fMgpPGR2WUJ4ncTERLKzs3t8PD4+nszMTIcdrzeZmZlkZ2eTm5s74NcO5bzDlQTUIJyub2bbwQqumj4WgNHB/lwzaxyvfHqUhma5KFIIV1q7di0ZGRlOP09+fj42m4309HSnn0sYJKAG4b0vymlp0yydMbbjvrSkiZyubyZ/z3ETKxNi+ElOTsZqtTr9PDabjfDwcAAsFovTzydkJYlB2bTnBGEh/lw0KazjvkviIxg/OoicgsPcOHeCidUJcX6/fnUXu4+aO046c8IofrVslqk1DFZ4eHjHTt5JSUkmVzM8SAtqgFrbNO/sO8GV08bi63N2+qSvj+K2hIm890U5ZacbTKxQCO9z4MAB0tLSCAsLIzExEZvN1vFY97Edm81GYmIiSilSUlJIS0sjPj6eNWvW9Ot4PUlOTsZms1FVVdWlBZWRkYFSioyMjH4dp11VVRUZGRmEhYV1GUdLS0vr0mW5Zs0awsLCurw2LCyMoqKiHo+dnZ1NRkYGVVVV5ObmkpubS0ZGRq+vcUtaa4/7SkxM1GYpOHhKx2a+pl/97Mg5jxWXn9Gxma/pv276woTKhDDs3r3b7BIcKiEhQVssFl1ZWdlxOzU1tcvja9eu7XI7KytLa631qlWrdHJy8oCO15PKykqdnJysV61adc79nc/f28/Rvc7Ox0pNTdXp6ek6JydHW63WjvuTk5N1cnKyLiws1FprnZeXpy0WS4/nycvL05WVlTo9Pb3Lz56Tk9Ovn9MZ+vqdBAr0eT7rpQU1QBv3nMDPR7FkSuQ5j8VFhLIgLpzcQrkmSghHWr58eUerZcWKFVRVVfX43KKioo6JDCtWrOjolhvs8cBo7dhsNnJyclizZk2XllJ+fv6AJ07k5+dTVVVFVlZWx33r1q0jOzu7o6XWfg6bzUZaWhrr168HIC8vj+XLl/d4bKvVisVioaCgoMvxO4+heQoJqAHatPcE8+PCGR18/qVkUpMmUnyylsKSShdXJoT3io+P73K7t8VHk5OTO7r81q9fT3Jy8pCOB0a3W0JCAhaLhVWrVnV0wRUVFZ33+H0pKio6Z2JHe2AWFBSQnJxMbm4u+fn5JCcnk5SU1DG1PTc3l7S0tB6P3X5cm81GQkJCx/15eXmkpKQMuFYzSUANwOHKOvaW1XSZvdfdDXPGExLgS27hYRdWJoToLC8vj/j4eGw2G+vWrRvSsfLz88nPz++4nZWVRUVFBZmZmVRUVAxqRp/Vaj2nZdfeiktKSiItLY28vLyOUElISKCiooKioiJsNlufoZifn3/ORI72sPMkElAD8I599Yj265/OJzTQj+tmj+e1z49R19TiqtKEEHY2m421a9dSWFhITk7OkKeEt7dIOk8wyMrKYs2aNV26Btu77fojNTWV8PDwjpZYVVUVaWlppKamYrFYWL58OQUFBeTn55OamgoY3ZKrV6/uuN2b7q2loqKijhagJ5GAGoCNe08wOSIUa+SIXp+XljSRM40tvLWzzEWVCSHaWSwW4uPjCQsLQylFWFjYgFaa6M5qtZKTk0NmZiZr1qwhOzubqqoqKisrycvLIyMjg8zMzHNm9/WlsLCwY++uyZMnk5CQQE5OTsfPYLVau4wZpaSkkJuby4oVK/o8ts1m6xJkPXV1ur3zzZww8wvI6us5Zsziq21s1lN+8Yb+zau7+nxuW1ubXpK1Sd++9kMXVCZEV942i28g1q5dqxMSEjpm6GmtdWFhoQY6ZsENRwkJCTovL8+083vFLD6lVDLg/EvCB2Hrl6doamljaS/de+2UUtyaEM2HtlMcr5ZrooRwte6THiwWi8d1bznSYCdzmM1tVpJQSlmB/l/l5mKb9h5nZKAfSXH9m6Z53ezx/Dn/C/L3HOeui2OdXJ0QAuiY7p2WltYxTTspKYmcnByXLIfkjoqKivo1buWO3CagAKvWOr+nza2UUulAOkBMTIwr60JrzcY9J7hsaiQBfv1rdE6NGkFMeAgbdklACeFK6enpsqBrJ53HtjyNW3TxKaWStdb5vT1Ha52ttU7SWidFRp57kawz7TpazYmaxl5n73WnlOLqmVF8eOAUZxplNp8QQgyUWwQUUKGUSlZKpQJWpVRCn69woY17TqAUXDFtYMGYMjOKptY2Nu8rd1JlQgjhvdwioLTWRfYWVDhgMbmcc2zae5x5kyyMGRE4oNclxoYRHhrAht0y3VwIIQbKLQKqnb0bL15r7TZL7p6oaeCzw6dZOiNqwK/18/Xhqulj2bT3BM2tbU6oTgghvJdbBZQ72lZsrKl36QURg3r91TOjqGlo4WNb72t9CSGE6EoCqg9FpZUE+vkwc/yoQb1+yZRIgvx9yJNuPiGEGBAJqD4UlVYyJ3p0v6eXdxcc4MviCyLJ231ctuAQQogBkIDqRWNLK7uOVJMQG9b3k3tx9awojp5uYJfJ220L4clyc3M71thr33qiN9132u1+u7vOu9o6q6beZGZmkp2dPajj9PWzeSp3ulDX7ew6Wk1TaxsJMZYhHWfp9LH4KNiw+zizo0c7pjghhhGbzcbKlSspLi522pJFa9euHdBqE46sKT8/H5vNRlZW1pCDzptIC6oXRfZNBxNihtaCGjMikKTYcDbsknEoIQajfX8jZ66nl5ycPKCAcmRNnXe7Hc5rBnYnAdWL7aVVRFuCGTsqaMjHSpkZxd6yGg5V1DmgMiGENwkPD+/YwLD7RoPDmXTx9aKotJLEIY4/tUuZGcXDb+xhw+7j3LN4skOOKcSAvPlTKNthbg3j5sB1vx/QS9r3YQIICwsjKyuL9PR0EhMTWbFiBatWrQJgzZo1rF+/nsLCwkGVlpiYSEZGRsc6fu3Hz8vLo6CgoGNfKKvV2mNNVVVVZGZm8vzzzxMeHk5qaipZWVl9njs5OZmVK1ees6dURkYG2dnZpKenk5mZ2e8WXk91pKWlER4eztq1awHjPVu9ejWVlZUdrw0LC2Pjxo1dtovvLDs7m8LCQrKysjp2Gm7fF6un1wyWtKB6cOx0PcdONwy5e69dXEQoU6NGyHRzIQYoKyuLtWvXkpCQQGVlpUsXgl29ejU5OTkdH+Dtkyh6qmnp0qVYLBYqKys5cOAANputY9fcviQlJbF69eou97WfZ6DjYz3VsWLFii7b1+fl5ZGUlNSxW3D7Yz0FTX5+PsuXLwfo2AE4NTWVlJSUc2p3BGlB9WB7aRXAkGfwdXb1zHE8+u6XVNY2ERYa4LDjCtEvA2y5CGOb9fYWTXtrqiftW753bjGtW7euo4XV09hSVVUVNpuNnJwcwsLCyMjI6Aij/Pz8AQdyX3XYbDZsNhtWqxWbzUZmZibr168nISGBvLy8jgA6H6vVisVioaCggHXr1nXc33kMzZGkBdWDopKhXaB7Pikzo2jTsGnvCYcdUwjhPPHx8V1ud98IsbOioqJzWjntodQ+vnQ+aWlpJCQkYLFYWLVqVUeLa7CbDPZVR3JyMrm5ueTn55OcnExSUlLHzMHc3FzS0tJ6PHb7cW02W5dWVl5eHikpKQOutS8SUD0Y6gW65zMnejTjRgWRt/u4w44phHAPVqv1nCCqqqoCep74kJ+f36XLLSsri4qKCjIzM6moqBjUjL6+6khLSyMvL68jVBISEqioqKCoqAibzdZnKLbPXux+nzN27JWAOo/GllZ2Hqlm3hCvf+rOx0eRPHMsm/eX09Dc6tBjCzHcWK1WTp06BRgfwOvXrze1ntTUVMLDwztaQFVVVR3jND0FTXuLpH0MCIyQWrNmTUeowNluO0fUsXz5cgoKCsjPz+/YaXf58uWsXr26Xzvvdm8tFRUVdbQAHU0C6jx2d1yg67jxp3YpM8dR39zK1i9POvzYQgwn7TPc4uPjSUtLc4vp2YWFhVRUVBAWFsbkyZP73M22fWZg+6zA7OxsqqqqqKys7JgZl5mZec7svqHUYbFYsFqtXcaMUlJSyM3NZcWKFX0e22azdQmy9evXO6X1BKA8cX24pKQk3Vuf7lD966MSHnx5J+9nXsnEsBCHHruppY3Eh/K4fs54slLnOvTYQgDs2bOHGTNmmF2GGCYSExPJysrqNaT6+p1UShVqrc/5C0NaUOexr6yGkYF+RFuCHX7sAD8fLp8Wyca9x2lt87w/DoQQorPBTuboDwmo89hXVsPUcSNRSjnl+FfPGsfJM018eqiy7ycLIYSbKioq6te41WBJQHWjtWZvWTXTxo102jmumBaJv69iwy6ZzSeE8Fx9jbENlQRUN2XVDVQ3tDDdiQE1KsifhdYxbJA9ooQQokcSUN3sLasBYFqU8wIKjK3gi0/WcqD8jFPPI4QQnkoCqpt99oCaPs5xK0icT/LMKADelm4+4QTSMhfuYii/ixJQ3ewrq2HcqCBGh/g79TzjRwdz0SQLr3x6RD5MhEMFBARQX19vdhlCAFBfX09gYOCgXisB1c3eshqnTpDo7Pb5k9h//AyFJTKbTzhOREQEhw8fpqKigubmZvkDSLic1prm5mYqKio4fPgwY8aMGdRx3GI1c6WUBbDav+ZrrTPNqKO5tY0DJ86wZEqES8637MIJ/Pb1PTzzcSlJcY5fCVgMT6NHjyYwMJDy8nJOnTpFS0uL2SWJYcjPz4+goCBiYmIIChrcpq9uEVDAcgCtdbZSar5SKl1rne3qIg6erKWptc3pEyTahQb6cfO8CTxfcJhfLpuJJUS24BCOERQUxKRJk8wuQ4ghcYsuPq11dqdAsgK27s9RSqUrpQqUUgXl5eVOqaNjBp+LuvgA7lwQS1NLGy8WHXHZOYUQwhO4RUC1U0pZgQqtdX73x+whlqS1ToqMjHTK+fcfr8HXR3HB2BH9e4HW8Nl6qB38wq8zJ4ziwkkWnvmkVMYKhBCiE7cKKCBVa92//ZGdwFZey6SwYIL8ffv3gs+eg5fS4aNHh3TeuxbE8OWJMxTIZAkhhOjgNgGllErVWq+xf5/Q1/OdoaSilpgxof178pkT8NZPje9tm4d03hsvHM/IQD+e+bh0SMcRQghv4hYBpZRKBrKUUoVKqULA5VPatNaUnKojNryf22u88RNoroNZt8DRImg4PehzhwT4cfO8aF7fcYzK2qZBH0cIIbyJWwSU1jpfax2vtU60f50zBuVsVXXN1DS0EDumHwG15zXY/TJcnglJ94Bug4Nbh3T+Oy+OMSZLbJfJEkIIAW4SUO6gtKIOgJi+WlD1VfD6jyBqDlx6P0xaAH7BUDy0br4Z40cxL8bCMx+XyGQJIYRAAqpDiT2gYvsag8r7T6g9AV/5X/D1B79AiFk45HEogDsWxHCgvJZPiiuGfCwhhPB0ElB2padqAZgU3ssuurbNUPQPWPR9mDDv7P3Wy6F8D9QMbeHXZXMnMDLIj2c/kckSQgghAWVXcqqOyJGBhAT0sLhGUx28+gMIt8IVP+v62OTLjf8WbxlSDcEBvtw6L5o3dpbJZAkhxLAnAWVXUtHHDL53HobKg7DsLxDQ7XnjL4Sg0VD87pDruMM+WeKFosNDPpYQQngyCSi70lN1xPQ0g+9IoXExbuI3YPKScx/38YW4JWDbYqwuMQTTx40iIUZWlhBCCAkooKG5lbLqBmLDzzNBoqUJXrkPRkRBym96Poj1CjhdCpXFQ67nzotjsZXX8rFMlhBCDGMSUMChjhl852lBbf0znNgFN/7J6MbrSfs4lANm8904dzyjgmRlCSHE8CYBBRyqNAJqUvcxqBN7YcsfYNatMO263g8SMQVGjh/y9VAAQf6+3Jowkbd2llEhkyWEEMOUBBRwuNLYHrvLFPO2Vvj3fRAQCtet6fsgShmtqOIt0NY25JruvDiGptY2XiiUyRJCiOFJAgojoAL9fIgcEXj2zk/WweFP4Nrfw4h+bu9hvRzqThldgkM0NWokSbFhPCuTJYQQw5QEFHC4so7osGCUUsYdlSWw8TdwQTLMXdH/AzlwHAqMlSVsJ2v5yCaTJYQQw48EFEYLKtpi797TGl57wOiyu/FPxn/7a3Q0jLnAIeNQADfMHc/oYH+ekZUlhBDDkAQUcKSynolh9gkSnz0LBzbB0l+BJWbgB5t8OZR8AK3NQ67LmCwRzVs7j3HqTOOQjyeEEJ5k2AdUXVMLp2qbmBgWbN+E8GcwaSHMv3dwB7ReDk1njIt7HeDOBTE0t2pyZbKEEGKYGfYBdcQ+g29iWPDZTQhv+l/wGeRbE7cEUA4bh5oSNZL5cTJZQggx/Az7gGqfYj6r+j37JoSrIHLq4A8YEg7j5zpsHAqMKecHT9Xx4YFTDjumEEK4OwmoyjpGUcvkj38JUbPh0geGftDJl8OhT6CpdujHAq6bbUyWeFpWlhBCDCMSUJX1PBjwLD515UbXnq//0A9qvRzamqH0w6EfC2OyxF0Xx/D6jmMUllQ65JhCCOHuhn1AhRzZynKfTahF34foBMccNGYR+Pg7bBwK4HtXXsC4UUE8+PJOWlqHvlKFEEK4u+EdUE11LD/2B8r8Jpy7CeFQBITCpAUOHYcKDfTjV8tmsudYNU99WOKw4wohhLsa3gH1zsOMbyvj3zE/PXcTwqGafDkc+xzqHLcKxLWzx3HFtEj+34Z9lJ1ucNhxhRDCHblNQCmlUpVSyUqpdJec8Egh+qNHeablKponXer4409eAmiHjUMBKKX4zU2zaWnTPPT6bocdVwgh3JFbBJRSKhVAa51vv53s1BPaNyFsDYlkdcudxjVQjhadCH5BcHCrQw8bMyaE7115Aa9/fowt+8sdemwhhHAnfmYXYDcfWG//3gYkAPk9PnvfPrjiikGf7KPZfiyM2MX/q7qbmqAQJj34YzhzbNDH69GFfpD3OPz+A4ceNkP58tLcb/DLR97mrc+fIEi3OvT4QgjRl4NBo6i98AiqNJKZx5wz5OAWLSjA0u32mO5PUEqlK6UKlFIFzc1DWedO4xvYyDuNCWxoS2JB9SGm1Z0cwvF6UTUCRtSDX4tDDxuoW3moOJ+DwWGsnbDAoccWQoje1Pr4s2bSEt5L9GdW8JecDgxw2rmUOyyfo5TKAvK01vn27r0UrXVmT89PSkrSBQUFQztpa7NjrnnqzcGt8OT1cPuzMP16hx/+vme38/auMjY8cBlxEaEOP74QQrTTWvPKp0dZ/eYeLj7zDn8J+Ct1CRmE3NSPDV37oJQq1Fondb/fXVpQ2zjbirICeU4/o7PDCTqNQ73vlMM/eMMMAnx9+OW/d8k6fUIIp/n0UBW3/d8HPLD+UxYEH+XPwX+HmEWE3PCwU8/rFgGltc4FrPbWk6V9soTH8w+CifPh4HtOOXzUqCB+dPVUtuwv582dZU45hxBi+Dp2up4frv+Umx/ZSmlFPX9eNom/qD/gExIOaU85/Q99d5kkgda6vZ3oHeHULm4JvLsa6ishOMzhh//awlhyCg7zm1d3c9nUSEYEus3/UiGEh6ptbCF7i421Ww7QpuE7V8TzvctiGLE+DWqOw7fehJFRTq/DLVpQXi1uMcb1UB855fB+vj48fMtsjtc08Oe8/U45hxBieGht0zz3SSlX/PFd/mfjFyydEcXGH15O5jXTGLHxZ1DyPnzlr8bwhQvIn9vOFp0IvoHGONS065xyinkxYdw+P4YnPjjIbYkTmTF+lFPOI4TwTlpr3tl3gt+/uZf9x8+QGBvG376aSGKsvdfnw0eg8ElY/EOYu9xldUkLytn8g4x1+Zw0DtUu89ppjA7258GXd9LWJhMmhBD9s720khXZH/GtJwtobGnj/+5KIPfbi86G07434e1fwIyb4Kr/dGltElCuELfYWJevvsppp7CEBPCz66ZTWFJJTuEhp51HCOEdDpSf4Tv/KuSWRz/AVn6Gh74yi/wfXs51c8ajlDKedKQIcr8FEy6CW9YOfqfxQZIuPlfoGIf60GndfACpiRPJKTjM6jf3kjJzHOGhzruATgjhmY5U1fM/+fvJLTxMkL8vDyRP4d4l1nMnWFWVwjMrICQC7ljv+AW1+0FaUK4QnXR2HMqJlFI8dPNszjS0kPXmXqeeSwjhWcprGvmvf+/iyj+8y8vbj/L1S+LY/JMreSB56rnhVF8JT6dBSyPcleOSGXvnIy0oV3DROBTAtHEjuWfxZNZusbF8/kQSY8Odfk4hhPs6XddM9nsHePz9gzS1tpGWOJH7lk4h2tLDItnNDfDsHXDqAHztRRg73bUFdyIB5Spxi2FzljEOFWxx6ql+sHQKr352lF+8tJPX7luMn680lIUYbk7XN/PY+8U88X4xNY0tLLtwAv+RPAVr5IieX9TWCi/eawxH3PYYTL7MdQWfhwSUq8ReCrrNuB5q2rVOPVVooB+/XDaLb/+rkCc/OMi9S6xOPZ8Qwn1UNzTz+PvFPPZ+MTUNLVw7axw/WDqFmRP6uPxEa3gzE/a8Ctf8DuakuqbgXkhAucrE+fZxqPecHlAA18yK4sppkfwpbz83zB3P+NFO2PNKCOE2qhuaeeL9gzz2vo3qhhaumRXFD5ZOYdaE0f07wHt/hG3rYNH3YdH3nFtsP/U7oJRSVwEVWutPO903GZgHFGmtDzq8Om/SsS6fcydKtFNK8eubZpPyp8385tXdPHpXwtmpo0IIr1Hd0MxTWw+y7j0jmFJmRnH/0inMju5nMAFseww2/RbmroCUh5xX7AD1GVD2YMoHtHFTaWCV1vq/tdbFSqlK4BTg69xSvUDcYtiyxiXjUGDsvnt/8hTWvLWPP27Yx0+uMW+wUwjhWCfPNPLE1mL+8UEJNY0tJM+I4oHkAQYTwM4X4fUfwZRr4CuPuPxap970GlD2FlI2kIERUhYgHliulPoCSNdav6PkT/P+iVsMm3/vknGodt+5PJ5DFfU88s4Bgv19+f5VU1xyXiGEcxytqid7i43ntpXS2NLG9bPH850r4gceTABfboQX0yFmIaQ96ZptiAagrxZUOpCotT7d6b7tQC6AUuo2pdRKvG0Fcmdx8TgUGF19D988m8bmVv64YT/BAX7cs3iyS84thHCc4pO1/O3dA7y4/TBaw83zovnOFfHE9zYrrzeHtsH6r0LkdLjjOVMuxO1LXwFV0C2cutBav2D/dp3jSvJi/kEwMQlKtrr0tD4+ijWpc2loaeWh13YT5O/DXRfHurQGIcTg7D5azaPvfskbO47h7+vDnQtiWHmZlYlhQwiUY5/D07fBiCj46gsuGXIYjL4CSlYddbS4xbDlD9BwGoIG0SQfJD9fH/68Yh4NzYU8+PJOgvx8uS1xosvOL4ToP6012w5W8rfNB9i09wQjAv1IvyyeexZPJnJk4NAOXr4P/nkzBIyEr//btFUi+qOvgJJlCByt/YLd0o9g6jUuPXWAnw+P3pXAPU9t4ye5nxHk78sNc8e7tAYhRM+aW9t4Y8cxHnu/mM8PnyYsxJ8fpUzl7kVxjA5xwPjQqQPw1E3g42eEkyVm6Md0or4C6gKl1IVa68+6P2Cf3ZcC2IC/aa1lFl9/TJwPvgHGOJSLAwogyN+XdXcn8fXHP+H+57YT5O/D0hnu+xeUEMNBdUMzz31SypNbD3L0dAPWiFB+e/NsbkuYSHCAgz5aqw7BP74CrU3wzTdgTLxjjutESuuee/GUUqOBYmAD3Wbx2W+v1FpXK6XatNYum5uYlJSkCwoKXHU6x3viemiug/R3TSuhpqGZr/79Y/Ycq+GxbySxZEqkabUIMVwdqqjjia0HWb+tlNqmVhZaw7l3sZWrpo/Fx8eBk6NryuCJ66D2lNFymnCR447tAEqpQq11Uvf7e21Baa1PK6WSgByMUALIA5K11tvtB74NqHJsuV7OpHGozkYG+fPUtxZwe/ZHrPxHAf/41sUsmCw9ukK4QlFpJY+9V8ybO4/hoxTLLpzAPYsnD26qeF9qT8E/boaa43D3y24XTr3ptQXV54uNFha9zfRzBo9vQRVvgaeWwZ3Pm9LN19nJM42sWPshx6sb+de9F3PRJIup9QjhrVrbNG/vKuPv79koKq1iVJAfd14cy9cviXXeUmR1FcaY06kvjG0zTF78tSeDakH1pT2YlFKjtNbVQznWsNIxDvW+6QEVMSKQp+9dyPK1H3L3Yx/zbPrC/q/dJYToU2VtEzmFh/jnRyUcqqgnJjyE/1o2k7SkSYR234fJkdrD6eR+uONZtw2n3gzp3VFKjQLWAPcO9VjDin+wsYmhi9bl68u40UE8fe/FLF/7IV977BOez1jIBWNHml2WEB5La01RaRVPf1TCazuO0dTSxoK4cH5x/UxSZkbh68jxpfPpHk4XLHXu+ZxkUKFiD6afAaswJlEM6d1WSlkAq/1rvtY6cyjH8whxi43VgxuqIaiPZfBdYFJ4CM+sNFpSd677mOczFhEXEWp2WUJ4lDONLby8/QhPf1zKnmPVjAj0Y0XSJO5aGMP0cS76d+4l4QQD3PJdKTVKKbUaqATSMCZOnNNvOAjLgSStdfsSSukOOKZ7i1t8dn8oNzE5IpSn772Y5tY27vr7xxypqje7JCE8wt6yah58eQcLf7eRB1/eiQJ+d8scPv75Uh66ebaE0yD1qwXVqcWUCRwAlrcvc9Q+UWIotNbZnW5aMWYKereOcagtMPVqs6vpMDVqJP+852LuWPcRd637iOczFjF2VJDZZQnhdhpbWnlzRxn/+qiEgpJKAvx8uHHueL66MJZ5kyyu397Gy8IJ+r4OahTwc4yuvELgp1rrjd2eMxpjn6ghX02mlLICmVrrjPM8lo6xeC0xMTGJJSUlQz2d+Z5aBtXH4PvbwM0WhC8sqeRrj31MtCWY59IXMmbEEJdXEcJLlJ6q4+lPSsgpOExFbRNxY0K46+JYUhMnEhYaYE5RHh5OPc3i6yugCjDW4zsnmDo9p18BpZRK5dylk2xa6/xOz1mltV7T23HAC6aZt9v2GLz+Q/jOBxA1y+xqzvHhgVN844lPiI8cwbMrFzpmqRUhPFBLaxvv7CvnXx+VsHl/Ob4+iuQZY/nqwlgujY9w7EW1A1V7ylghwkPDCQYfUKOBZKBSa72pl+cMuQWllErtNAaVoLUu6um5XhNQZ8rhv6fCkh/BVQ+aXc15bd5fzsqnCpgxYRT/vGcBo4IkpMTwsftoNS8WHeaVz45SXtNI1KhA7lgQw+3zYxg32g26vmvKjHCqPAi3P+OR4QSDDKhOL+4xqBwRUEqpZGAtZ1ekyOzcsurOawIK3Lqbr92GXWV85+kiokYG8uuvzCZlpqzdJ7zX8eoGXt5+hJe2H2FvWQ3+voorp43ltsSJLJ0+Fj9fN9lxtqrU6NarLTf2c5q8xOyKBm1IAdXpIKMxZtwdaA8qR45B9ZdXBVR7N9+3t8K42WZX06Oi0kp+9sIO9h2v4dpZ4/ivm2a5x1+QQjhAbWMLb+8q46XtR9j65UnaNMyLsXDrvGhunDvBvLGlnpz8Ev5xEzSdga++aOwz58EcElCdDtYeVKeAjUhADZ4HdPO1a25tY917Nv4n/wv8fX34yTXT+OrCWOdfdCiEE7S2aT44cJKXio7w1q4y6ppamRQezC0XRXNLwkQmu+t1gGU7jf2ctDbW1hs3x+yKhsyhAdXpoKMxZtZlyWrmQ/DUMqg+Ct8vcNtuvs5KT9Xxi5d38N4XJ7lwkoXVt8xh5gTzLzYWoj/2llXzYtERXvn0CMerGxkZ5MeNcydwa0I0SbFhrp8ePhCHC+Fft4J/iLEqecQUsytyCKcEVKeDj3blgrFeF1AFj8Nr/+H23Xydaa3592dHeei13VTWNXPv4sncnzyFkABZ8Uq4nxPVDbzy6VFe3H6EPceq8fNRXDFtLLcmRHPV9LEE+XvAdnYH34dnVkBoBNz9CoTFmV2Rwzhlsdh2rl7N3OtMXwav/wh2veQxAaWU4isXRXP51Eh+/+Ze1m6x8fqOYzx082yunDbW7PKEoK6phQ27jvPi9iO8/0U5bRounGTh1zfN4sa54z3r2r4v8mH9XWCJNcJp1PDYCdshLShX87oWFBizcaqPeEw3X3efFFfwsxc/50B5LTfOHc8vl81k7EiZRCFcq7K2iU17T5C3+zhbviinrqmVaEswt8yL5paEaOIjR5hd4sDtfgVy74GxM+BrLxktKC/j1BaUcIBZNxvdfMd3euSg54LJ4bxx/xLWbrbx101fsmV/OT+9bga3z59k7kWMwuuVnqpjw+4y8nYfp6CkktY2TdSoQG6ZF82yCyewIC7cc38HP30WXvmusTTanc9DsMXsilxKWlDuovYk/HEKLP4hLP1Ps6sZElv5GX7x0k4+tJ0iKTaM3906h6lRsn2HcIy2Ns2OI6fJ232cvN3H2Xe8BoDp40aSMjOKlJlRzIke7d6THfrj47Xw5ipjH6fbn4VAD2z99ZNTJ0m4mlcGFBjdfKcPw32FHtnN15nWmheKjvDw67upaWgh43Ir9101xTMGo4XbaWxp5YMDp8jffZz8Pcc5Xt2Ir49iflwYKTPHkTIjipgxIWaX6RhaQ/5/wdY/w7QbIPVx8Pfu7nLp4vMEs26B1x7w2G6+zpRSpCZO5KrpY3n49T088s4BXv/8GL+9eQ6Lp3hfH7pwvNN1zWzaZ7SSNu8rp7aplZAAXy6fGknKzCiunDbW/S6gHaqWJvj3ffD5c5D4Tbj+j+A7fD+mpQXlTmpPwh+nwuIHYOkvza7GoT748iS/eHknxSdruWVeNA/eMMOzZlEJlzhUUdfRdffJwQpa2zSRIwNJnhHF1TOjWBQ/xntb4Y018PzdcGATXPkgXPZjj+9J6S/p4vMU//gKVB3yim6+7hqaW3n0nS/5v80HCA304+fXzeCWhGj83WVtM+FyWmt2Hqkmb3cZG3YfZ2+ZMZ40ZeyIjvGkCydaPHeSQ3/VHIdn0oxVIpb9DyR8zeyKXEoCylMUPGF082W8B+Pnml2NU3xxvIafv7SDbQcrsYT4c83McVw/dzyXxI+RsPJyrW2aPceq2Xawgm0HK/ikuJKTZxrxUZAUG07KzCiSZ0a57zJDznDyS2N1iNpySHsSpl5jdkUuJwHlKby4m6+ztjbNO/tO8Nrnx8jbfZwzjS0dYXXD3PEskrDyCg3NrXx++LQ9jCooKqmkprEFgGhLMPPjwrj0ggiumj52eHb5Hi6AZ5Yb39+ZAxMTza3HJBJQnuQfN0NVCdxX5HXdfOfT0NzKlv3lvLHDCKvaplYsIf5cO2sc18+RsPIk1Q3NFJZUsq3YaCF9dvg0TS1tgNFtN39yOAviwpk/OZxoS7DJ1Zps/9uQ8w0IjTQuwB0Tb3ZFppFZfJ5k1s3w6v1QtsNru/k6C/L35epZ47h61riOsHp9xzFe/ewoz207RFiIP9fOtoeVdYz77McjKK9p7GgdbTtYwZ5j1bRp8PVRzI4ezdcXxTI/LpykuHDCvW3G3VAUPmVcmD9uNtyVCyNkebDzkRaUO6o9ZVy0e+n9kPwrs6sxTUNzK5vtLat8e8sqPDSAa2ZFccOcCSy0hktYuZDWmtKKuo4w2nawkuKTtQAE+fuQEBPG/LhwFkwOZ16MRRYOPh+tYfMaePd3EH8VLP8HBMpF7NLF52mGWTdfXxqaW3l3nz2s9hynriOsxnHj3PFcPFnCytFa2zT7ymqMFtLBCrYVV3CiphEAS4g/SbHhLJhshNLs6NHSDduX1hZ440dQ+CTMvR1u+l/wk1YlSBef55l1C7z6Ayj7HMZfaHY1pgvy9+Xa2eO4dva4jrB6fccxXvn0CM9+UsqY0ACumT2OG+ZIWA3UmcYWDp6sxXayluLyWopPnqH4ZC0Hyms5Y5/QMGF0EIvix3S0kC6IHOH9U78dqakOXrgH9r1hX87sl/KHZz9IC8pdSTdfvxhhdYLXd5Sx0d6yag+rq2dGYY0YwXhL0LD/676ppY1DlXX2ALKHkT2Ijlc3djxPKZgwOhhrZCiTI0KZF2Nhflw4E8O8ZBkhM1QfhWfvgGOfwXVr4OJ0sytyO9LF54n+eQtUFMMPtstfW/1Q39QeVsfYuOcE9c2tAPgoiBoVxMSwYKItwUwMCyE6LLjj9gRLsFesTtDWpjle00BxeS0HurWGDlXW09p29t96eGgAkyNCO76sEaFYI0cQOybEK94Lt3G4EJ67E5rOwK3rYPr1ZlfklqSLzxPNvFm6+QYgOMCX6+aM57o546lvamV7aSWHK+s5XFXP4co6jlTWs+1gJa9+fqzLhzVA5MjA8wbYREsw0WHBbjXgX1XX1Kk77myL6ODJ2o5QBmPiwuSIEcyaMJob507oaBVNjgjFEiJjH063Ixde+Z4xQ+9rGyBqltkVeRz3+VcnzjVjmTEVdddLElADFBzgyyUXnH9R2pbWNsqqGzhSWc/hynqOtAdYVT07jpzm7V1lNLd2DbDw0IBOAdY1yKJGBdHapmlsaaWhuY3GllYaW9pobG6joaWVxi732f/b0kZDc9f7Gjo91uVYzWefX9/U2nGhKxjTuSeFBTM5IpRL4sd0tIYmR4YSNTJIxonM0NYG7zwM7/0RYi6BFf/0yk0GXUECyp2FhIP1ctj1Miz9lXTzOYifrw8Tw0KYGBbCxed5vK1Nc6KmkSNVdUYLrFOQ7Ttew6a9J2i0X3w6FEpBoJ8PgX6+BPkb/w308yHQ34cg+32jg/3tz/EhyN+XIH9foi1GIE2ODGVSWAgBfsN7fM2tNJ6BlzJg72uQcDdc/98yU28I3C6glFJZWutMs+twG7NuMZbfP/YZTLjI7GqGBR8fxbjRQYwbHURi7LmPa605eaapo+V1oroRf3uItAdJe/C0h02gv8859/n7Ks/fVE+cVVVqTIY4sRuu/T1c/G35o3KI3CqglFLJgNXsOtzK9Bvh1Qdg98sSUG5CKUXkyEAiRwZy0SSL2eUId1D6ETx3F7Q2w105cEGy2RV5BbfpG1BKWQFbL4+nK6UKlFIF5eXlLqzMZCHhYL3CGIfywBmXQni97U/DkzdC0Ci4N1/CyYHcJqAAq9a6x4DSWmdrrZO01kmRkZGurMt8s26GyoNGN58Qwj20tcLbv4BXvguxl8C9GyFyqtlVeRWXdfEppVKB8G5327TW+UqpZK11vqtq8TjTbzw7m0+6+YQwX8NpeOFe+GIDLEiHa34Hvv5mV+V1XBZQWuvcXh6usI8/WQCrUipBa13kmso8QHs33+fPw5U/B79huG+OEO6iwgbP3A4VB+CG/wfz7zG7Iq/lFl18WusiewsqHCOkRHeX3Ac1R40dd4UQ5ijeAuuugtoTxh5OEk5O5RYB1c4+zhQvrafzsF4BcUuMi/+aas2uRojhp+BxY/mxEVGwchNMvszsiryeWwWU6MPSX0JtOXy81uxKhBg+mhuMSz1e+w9jD6d78iBcroZxBQkoTzJpAUy5Brb+D9RXmV2NEN7v1AF4LBkKn4BLH4A7njOmkwuXkIDyNFc9CA1V8OEjZlcihHfb9RKsvRxOH4Y71kPKr8FHVnp3JQkoTzN+rrHK+UePQu1Js6sRwvu0NMLrP4acb8DYGZDxHky71uyqhiUJKE905S+guQ7e/5PZlQjhXSps8FgKbFsHi74P33wDLJPMrmrYkoDyRJFTYe7t8Mk6Y7dOIcTQ7X7F6NKrPAi3PwvXPCwX35pMAspTXZEJug22/MHsSoTwbC2N8GYmPH83REwxuvRk51u3IAHlqcLiIPHrUPQPY1t4IcTAVR6Ex6+Fj/8GC78L33wLws6zx4owhQSUJ1vyY/Dxg81ZZlcihOfZ8xqsvcyYSr7iX3Dtatlc0M1IQHmyUeNhwUr4fD2U7zO7GiE8Q0sTvPVzWH+XccHtt7fAjGVmVyXOQwLK0136H+AfCu88bHYlQri/qlJ44jr46BFYkAHfetvoLhduSQLK04WOgUXfNWYgHf3U7GqEcF/73oS/LYGT+yHtKbh+jewM4OYkoLzBou9BkAU2/dbsSoRwP63NsOFBePZ2YwJExmZjE1Dh9iSgvEHQaFj8AHyZB6UfmV2NEO7j9GF44nr44H9h/r3wrQ2y0KsHkYDyFgvSjW0ANj4EWptdjRDm278B/rYYTuyB1Mfhhv8G/yCzqxIDIAHlLQJCjWnnJe+D7R2zqxHCPA3VxvYYz6TBqIlGl97s28yuSgyCBJQ3Sfw6jJ4krSgxfH25ER5dBIVPGrtQ35sHY+LNrkoMkgSUN/ELhMsz4WgR7HvD7GqEcJ2G0/Dv++Bft4J/MNyzAa7+rfG98FgSUN7mwjtgzAWw6WFoazO7GiGc74t8o9W0/V9w6f3w7feMzT2Fx5OA8ja+fnDFz+DELtj1otnVCOE89VXwyvfg6dsgYISxFXvKb6TV5EUkoLzRrFsharaxukRrs9nVCOF4+zcYraZPn4HF/wEZW2BiktlVCQeTgPJGPj7GpoYVNuMfsBDeor4SXvqOMUMvaDTcmw/J/yXTx72Un9kFtFNKJQBWAK11rsnleL5p10F0EmxeAxfeLku6CM+37y147QE4c8K4pOLyVfJ77eXcqQX1M3swhSul5FLvoVIKlv4nVB+GgifMrkaIwaurgBcz4NkVEBwGKzcav9sSTl7PLVpQSql0YJtSyqq1zja7Hq9hvQLilhi77s68CUZNMLsiIQZm7xtGq6n2JFy2Ci77sQTTMOIuLah4YAxQoZRaq5SydH+CUipdKVWglCooLy93eYEe67o10NIAz6yAxjNmVyNE/9RVwAsr4bk7IDQSVm6Cq34h4TTMuCyglFKp9pDp/JXc6SkHtNZVQCGQ3v31WutsrXWS1jopMjLSVWV7vqiZkPoEHN8JL9wLba1mVyRE7/a8Bo9cbFwmcflPYeU7MOEis6sSJnBZF18fEx+2AeH27y1AlbPrGVamXm20pN74sbHtwLWrza5IiHPVnoI3V8HOXIiaA199AcbPNbsqYSK3GIPSWucqpVa1t6hkHMoJFqyEUwfgo0eN7QYWrDS7IiHO2v0KvP4j4+LbK34OS34Ivv5mVyVM5hYBBaC1XmP/Nt/UQrzZNQ9D5UHjr9SwOJiSYnZFYrirPWm07He9BOMvhK+9DONmm12VcBPuMklCuIKPL9z2d2OViZxvQNlOsysSw1VLI3z4CPxvojHmdNWDcO9GCSfRhQTUcBM4Au5cD4GjjJl9NWVmVySGE61h18vwyAJ4++cwYZ6xuOtlP5EuPXEOCajhaNQEI6TqK42Qaqo1uyIxHBz6BB67GnK+Dn7BcNcLcPfLMHaG2ZUJNyUBNVyNn2tsg132ObyYLtPPhfNUFBtdyo+lQFUJLPsLfPt9mJLc50vF8CYBNZxNuxauWQ17X4O8X5pdjfA29ZXw9i+M7rx9bxmbad5XZOz87Os287OEG5PfkuFu4beh4gB8+Fdja+ykb5ldkfB0LU1Q8BhszjKmjV90l7EKhCy1JQZIAkoYraiKYnj9x2CJhQuWml2R8ERaw55/Q96voLLYWAvy6t/CuDlmVyY8lHTxCaO7Je0JY7A65xtwfLfZFQlPc7gAHr8Wnr/bWC/vrlz7NU0STmLwJKCEIXCkMbPPP8SY2XfmhNkVCU9QeRByvwV/X2p0Fd/4Z/j2VuMicKXMrk54OAkocdboiXDnc1B3Ep69HZrqzK5IuKv6KmNdx7/ON7bEuOwn8IPtkPRNmQAhHEYCSnQ1YZ6x2sSRIngpA9razK5IuJOWJvjob/CXi+CDv8KcNLiv0FgJInCk2dUJLyMBJc41/QZjcHvPv2Hjr82uRrgDrWHPq/DoQngr0xhbytgMNz8Ko6PNrk54KWmLi/Nb9D1jTGHrn43VzxO/bnZFwixHCuHtB6H0A4iYBnc+D1OuljEm4XQSUOL8lILr/gCVJfD6DyEs1pg2LIaPqlLY+BvYkQMhEXDD/4MEuchWuI78pome+fpB2pPw+DWw/m742kswMdHsqoSzHSmEj7ONHW2VDyz5MVx6PwSNMrsyMcxIQIneBY0yunQevwYeS4aF34Urfmasii68R0ujscr4J2uNgAoYYbSWFj9gzO4UwgQSUKJvlknwnQ+MCRMf/hV2/xtu+G9jK3nh2U4fgcInoPBJqC2HMVOMrt0Lb5cWkzCdBJTon2AL3PgnmLMcXr0fnkmD2bfBtb+HEWPNrk4MhNZQ8gF8km3MzNNtMPVauDgdJl8BPjK5V7gHCSgxMLGLjA3m3v8zvPdH+DLfmJI+72syq8vdNdXBjufhk3VwfCcEWWDRd2H+vRAWZ3Z1QpxDaa3NrmHAkpKSdEFBgdlliPL98NoDULIVYhfDsj9DxBSzqxLdVRTDtr/D9n9Cw2mImg0L0o2LbANCzK5OCJRShVrrpO73SwtKDF7kVPj6a8YHX95/wv9dYix5c+kD4BdgdnXDW1sb2N4xWkv73zJm4828yQimmEXS2hUeQQJKDI2Pj3ER79Rr4a2fwjsPw84XjF1TYy42u7rhp6EaPnvWCKZTX0BopPFHQ9I3ZT8m4XEkoIRjjIwytuy48HZ4/Ufw+NWQdA8k/wqCRptdnfcr329MevjsWWg6A9FJcOs6mPkVY/sLITyQ2wSUUioVqAKsWutsk8sRgzX1Goi9FN75HXz8f7D3dbj+DzBjmXQrOVpbK+x/27h2yfYu+AYYMysXrIRouaBaeD63mCShlEoGKrTWRZ2/7+n5MknCQxwpgld/AGU7YNoNRlDJwqJDV1dhjPtt+7uxHNGoaEj6FiR+A0IjzK5OiAFz90kSBUChUioNowWVb3ZBwgGiE2Dlu/DRo0aL6pEFsPRXMP8e8PE1uzrPcqYcDr4HX+QZSxC1NEDcEmOK/7QbZH084ZXcogUFoJRaBWQA+VrrjPM8ng6kA8TExCSWlJS4uEIxJJUH4bUfwoGNxvjITX+BqFlmV+W+6quM6fvFW6D4PTixy7g/YCTMuc2YjSfvn/ASPbWgXBZQ9jGm8G5327TW+d26+LKAbVrr3J6OJV18Hkpr2JFrzPZrqIJF34eL7oSIqTI+1VQLpR/aA2kLHPvMWOHBLwhiFsLky2Dy5TD+ImktCa9jekD1Rim1Smu9xv69BVje20QJCSgPV1dhbBf+6dPG7dBIiL3EuNg39hIYO9P7l9tpboDD284G0pECaGsBH3+YON8eSJfBxCSZhSe8nrsHlAVYDtjoxyw+CSgvUWGDg+/Dwa1Gd9bpQ8b9wWEQc4kRVnGXwri5nj9m1doMR7dD8Wajy+7Qx8Y4kvKBCfPOBtKkiyEg1OxqhXAptw6ogZKA8lJVpfawet9YzLTCZtwfOMro5mpvZU24CHz9TS21T21tcHzH2RZSyQfG9UkAUXPOBlLsIrlOTAx77j6LTwiwxMBFMXDRHcbt6qPGB3vJViO4vthg3O8fCpMWGNdbxV1qXPNjdjeY1lC+zx5Im42WYUOV8VjEVOMC5smXGQEbOsbUUoXwFNKCEp7jTLkRViVbjeA6vtO43y/IGLeJvcQIrYnzB74IqtbQXA/NdcZXU93Z75vrjUkMzfXQbP9v58fPnDDqqT1hHMsSe7aFFLcERo137PsghJeRLj7hfeoqjJlvJR8YLZayz42Zbz7+xjVY4y8yJh50CZ0eQqa5buDn9w00gjBoNExcYA+lJbJ1hRADJF18wvuEhMP0G4wvMLaSKP34bCvr02eMrr+AEPDv9BUaCf7BRlehf3DXxwPs93Xc7vza4K6Pe/rEDSHcnASU8B5Bo41t6GUreiG8gpdfbCKEEMJTSUAJIYRwSxJQQggh3JIElBBCCLckASWEEMItSUAJIYRwSxJQQggh3JIElBBCCLckASWEEMIteeRafEqpcmCoe75HACcdUI63kPfjLHkvupL3oyt5P7pyxPsRq7WO7H6nRwaUIyilCs63OOFwJe/HWfJedCXvR1fyfnTlzPdDuviEEEK4JQkoIYQQbmk4B1S22QW4GXk/zpL3oit5P7qS96Mrp70fw3YMSgghhHsbzi0oIYQQbmxYbFiolEoFqgCr1vqc5mhfj3uT3n5WpZQFsNq/5mutM11eoIv19/+9UipL3g9QSiVg/H6gtc51bXWuJ58dZ9l/1gytdUovj1fhwPfC61tQ9jcNrXW+/XbyQB73Jv34WZcDSe0fPEqpdNdW6Fr9/X9vv9/qwtJM0c/342f2349wpZRXvyf9+OxIBmz2x2328PZavf1B4qzPUa8PKGA+YLN/bwO6/xL19bg36fVn1Vpnd/rLx9rpud6qz//39g9hb38f2vX6ftj/YNmmlLLaf1e8/X3p6/ejAMhpb1VqrYtcWZybccrn6HAIKEu322MG+Lg3sXS7fd6f1f6hXNH+15AXs3S7fb73wzoMPojbWbrd7v5+xNvvq1BKrbV3CXszS7fbXd4PrXUVsBbIARJdU5LbsnS77ZDP0eEQUFVA+BAe9yZV9O9nTdVaZzi5FndQRS/vh1IqeRiEdGdV9P37ccD+wVwIeHUXMP34/QDytdbxQFV7N9cwVYUTPkeHQ0Bt42y6W4G8AT7uTfr8WZVSqVrrNfbvvbm7E/p+PyqUUsn2Dx6rvB9s6/S9BeNDyZv19X4kdOrWW83w+UP3fJzyOer1AWUf2LPa/9qxdBrEy+vtcW/U13thvz9LKVWolCrEy//B9eN3o8h+XzjndmF4nX7+W7G0D4B7+6y1vt4PIFsplW5/fLm3vx/2nzOpc0vR2Z+jcqGuEEIIt+T1LSghhBCeSQJKCCGEW5KAEkII4ZYkoIQQQrglCSghhBBuSQJKCCGEW5KAEkII4ZYkoIQQQrglCSghhBBuSQJKCCGEW5KAEsKNKKVSlVIHlFKV3b4OKKWyzK5PCFcaFlu+C+EJ2neo1VrHK6XStdbZ7SuoD/PN8MQwJS0oIdxHRadttS32/yYzfHb0FaILCSgh3IR9I8COlpRdfPv9Qgw3ElBCuJ9MoH0/HWtvTxTCm0lACeF+kjuPObVvECjEcCMBJYQbse9WmtvpriKGwW6+QpyP7KgrhBDCLUkLSgghhFuSgBJCCOGWJKCEEEK4JQkoIYQQbkkCSgghhFuSgBJCCOGWJKCEEEK4JQkoIYQQbun/A8kO6wB179/5AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# To make figure 6\n",
    "\n",
    "# linear interpolation \n",
    "index = next(i for i, x in enumerate(pi_array_3) if x > pi_array_0[1])\n",
    "print(index)\n",
    "\n",
    "Omega_array_3_new = np.zeros(len(pi_array_0)) #living on the grid of pi_array_0\n",
    "\n",
    "for j in range(1,len(pi_array_0)-1):\n",
    "     if pi_array_3[index-2+j] < pi_array_0[j] <=pi_array_3[index-1+j]:\n",
    "            Omega_array_3_new[j] = ( Omega_array_3[index-2+j]+(Omega_array_3[index-1+j] - Omega_array_3[index-2+j])\n",
    "                         *(pi_array_0[j]-pi_array_3[index-2+j])/(pi_array_3[index-1+j]-pi_array_3[index-2+j]))\n",
    "\n",
    "Omega_array_3_new[0] = np.nan             \n",
    "Omega_array_3_new[-1] = np.nan\n",
    "pi_array_0_new =list(pi_array_3[:index-1])+list(pi_array_0[1:])\n",
    "\n",
    "diff = Omega_array_3_new - Omega_array_0\n",
    "diff_new = list(Omega_array_3[:index-1]) + list(diff[1:])\n",
    "diff2 = aveOmega_0 - Omega_array_0\n",
    "\n",
    "plt.axhline(y=0, color='r', linestyle='-')\n",
    "plt.plot(pi_array_0_new,diff_new,label=r'high $\\&$ low $\\eta$')\n",
    "plt.plot(pi_array_0,diff2,label=r'full info $\\&$ low $\\eta$ ')\n",
    "plt.xlabel(r'$\\pi$', fontsize=16)\n",
    "plt.ylabel('$\\Delta \\Omega$', fontsize=16)\n",
    "plt.legend(fontsize=16)\n",
    "\n",
    "plt.tight_layout()\n",
    "\n",
    "# DataOut = np.column_stack((pi_array_0_new,diff_new))\n",
    "# np.savetxt('usability_diff_highloweta_lowalpha.dat', DataOut)\n",
    "# DataOut = np.column_stack((pi_array_0,diff2))\n",
    "# np.savetxt('usability_diff_fullinfoloweta_lowalpha.dat', DataOut)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
